基于线性规划的建厂问题

基于线性规划的建厂问题
摘要
随着我国经济的不断发展，特别是工业的现代化，关于如何建立工厂并使其获
得最大的经济效益这一类问题成为工业经济方面的常见问题。

本文研究的是建厂方案问题。

内容为确定一个建厂的方案（包括建厂的地点，生产规模，运输销售方案等各因素），使得超轻型合金销售的年平均利润
最大，其实质就是利用线性规划求解最优解问题。

针对第一问求三市之间距离，本文建立了几何数学模型，对三市简化为地
球表面上的点，地球看做标准球体使得问题简化。

利用弧长公式以及MATLAB软件即可求得：A1、A2两市距离为132.6625 km; A1、A3两市距离为309.5062 km; A2、A3两市距离为294.8311 km。

针对第二问建厂方案问题，建立了线性规划模型。

将整个建厂方案简化为：年平均利润=销售额-原料成本-年均建厂成本-运输费用-人工费用-废料处理费-保证金。

问题综合考虑产品售价，原料成本，建厂成本，运输费用，人工费用，废料处理费，环保保证金，原料限制等因素，列出利润的目标函数以及相关约
束条件。

针对模型的求解，本文应用MATLAB软件中线性规划的相关知识求解，在得出目标函数和约束条件后，利用MATLAB软件编写目标函数与约束条件关系便可得到满足最优解的各变量的值以及最大年平均利润。

其建厂方案为：A1,A3作为建厂地点，A1生产规模为192860吨/年，A3生产规模为96430吨/年；对于主要原材料ST21运输，A1市的310000吨全部用于本地生产，A2市的230000吨全部运输到A1生产，A3市的270000吨也全部用于本地生产；对于超轻型合金的销售，A1市本地销售120000吨，A1市运输到A2市销售72860吨，A3市生产的96430吨全部运输到A2销售；这样可以达到最大年平均利润2.2184*10^9元。

在模型的评价与推广中，将所得的各项数据联系实际情况，考虑模型的优缺点，得到数据符合客观事实，具有一定的真实性和稳定性。

然后对模型进行了扩展，进而得出工业类问题的最大经济效益模型。

关键词：球面距离，线性规划，MATLAB软件
一、问题重述
生产某种新型金属合金‘超轻型合金’，主要原材料是矿石ST21，主要产
地是A1、A2、A3市，三地ST21产量分别为31万吨/年，23万吨/年，27万吨/年，它的运输费用为4.8元/吨*公里。

超轻型合金在生产时还需要耗费其他原
材料，以及不同地方不同的人工费和废料处理费（见附表1）。

A1、A2两市建
有超轻型合金市场，市场需求分别为12万吨/年，17万吨/年，超轻型合金的
运输费用为5.2元/吨*公里。

A1、A3市建厂，生产规模不受限制；如果在A2
市建厂，生产规模受到环保新规定的相应限制。

建厂成本和生产规模成正比，
比例系数为ki（见附表1)，使用年限50年。

问题一：根据附表1，计算A1、A2、A3之间的距离（以球面距离计算即可），并用MATLAB绘图工具绘制示意图。

问题二：设计建厂方案（建厂地点、生产规模、运输销售方案），使得超
轻型合金销售的年平均利润最大。

建厂地点的选择就是确定工厂坐落的具体地点，工厂建在什么地点，不仅
影响生产规模，而且还影响工厂的运输销售。

它直接关系到一个公司在未来的
销售利润或者得失和成败。

附表 1
二、问题分析
问题一：已知A1、A2、A3经纬度，计算三市之间的距离。

这属于求解球面上点之间的球面距离的几何数学问题，对于解决此类问题一般根据L = nπR/180，结合MATLAB中distance函数即可求解。

MATLAB绘制示意图：三维图真实反映了三市实际位置，但是A1、A2、A3
之间距离相对地球表面较小，无法清晰显示三市位置关系。

因此结合二维图，
以经纬度为纵横坐标可弥补上述缺点。

问题二：选择建厂位置，合理安排生产规模与运输销售方案，以取得“超
轻型合金”最大年平均利润。

它属于线性规划数学问题，对于解决此类问题，
我们必须首先确定目标函数，其次是约束条件。

经分析，此题目标函数：年平均利润=销售额-主要原材料ST21费用-其他
原材料费用-人工费-废料处理费-排污保证金-原材料ST21运输费-超轻型合金
运输费-年均建厂成本。

其约束条件受到多方面影响，如：三市ST21的产量；
A1、A2市场需求；A2生产规模等。

表示约束条件的时候，不能改变各变量范围，相应变量对应相应界限；此模型约束条件不仅包括不等式，也含有一系列等式
约束。

只有全面考虑，才能使得结果符合实际情况。

三、模型假设及符号说明
3.1 模型假设
①地球是标准球体；
②A1、A2、A3市相对地球表面（参考物）均看做点；
③建厂后50年内主要原材料ST21价格、其他原材料费用、人工费、废料处理费、ST21运输费、超轻型合金运输费均保持不变；
④建厂后50年内三地ST21产量、A1,A2两市市场需求保持不变；
⑤建厂一次性可确保使用50年。

3.2符号说明
四、模型的建立与求解
4.1问题一：计算A1,A2,A3之间的距离(球面距离)
A1、A2、A3三市经纬度：
地球半径为6371公里，根据几何关系球面上任意两点之间的球面距离：
L = nπR/180⑪式
利用MATLAB中distance函数可以得到的球面上两点圆心角度数
n = distance(纬度1，经度1，纬度2，经度2)
代入⑪式，结合MATLAB 语言编程(见附录1)得到：
A1、A2、A3三市之间距离：
图1：
A1、A2、A3三地位置关系三维图 (见附录2)
图2：
纬度
经度
A1、A2、A3三地位置关系二维图（见附录3）
4.2 问题二：建厂方案设计
4.2.1 模型的分析
根据“问题分析“的结果，需要利用线性规划的知识建立利润最大的数学模型，使得超轻型合金销售的年平均利润最大。

4.2.2 建立线性规划模型
4.2.2.1 目标函数的建立
由上述“问题分析”得到：年平均利润=销售额-主要原材料ST21费用-其他原材料费用-人工费-废料处理费-排污保证金-原材料ST21运输费-超轻型合金运输费-年均建厂成本，则有目标函数：
max z=(31515-3694.5-2.8*5000-908.654-3509.8)S1+(31515-3694.5-
2.8*5000-1439.6-442
3.7)S2+(3151-369
4.5-2.8*5000-1317.8-2721.6)S3-
4.8*132.6625*(Y12+Y21)-4.8*309.5062*（Y13+Y31）-4.8*294.8311*（Y23+Y32）-
5.2*132.6625*(Q12+Q21)-5.2*309.5062*Q31-5.2*294.8311Q32-4423.7*S2*50%-31996*S1/50-48187*S2/50-36425*S3/50
简化，得：
max z=8763.126S1+4781.61S2+9052.6S3-636.78Y12-1485.6297Y13-636.78Y12-1415.1892Y23-1485.6297Y31-1415.1892Y32-689.845Q12-689.845Q21-
1609.4322Q31-1533.1217Q32
4.2.2.2 约束条件的建立
由“问题分析”可知：
1. A1、A2、A3三地ST21产量分别为31万吨/年，23万吨/年，27万吨/年，本地拥有量与外地运来量之和大于生产用去量与运到外地量之和，可得①
②③。

2. A1、A2、A3三地生产规模受到原材料ST21总量与市场需求限制，而29万吨的市场需求大于原材料ST21总量可生产的超轻型合金，所以应以后者来控制对变量的影响，可得⑥。

3. A2建厂，规模不超过6万吨，可得⑦。

4. A1、A2市场需求分别为12万吨/年，17万吨/年，超轻型合金在这两地的销售应不超过市场需求，可得④⑤。

5. A1、A2、A3三市生产规模等于本地出售量与运到外地出售量，可得⑧
⑨⑩。

6. 各决策变量非负，可得11。

从而得到以下约束条件：
4.2.3 模型的求解
利用MATLAB求解线性规划（见附录4）,得到结果：x =
1.0e+005 *
1.9286
0.0000
0.9643
3.1000
0.0000
0.0000
2.3000
0.0000
0.0000
0.0000
0.0000
2.7000
1.2000
0.7286
0.0000
0.0000
0.0000
0.9643
fval =
-2.2184e+009
4.2.4 结果的分析
变量结果（未说明的全部等于零）：
由此可知，应选择A1,A3作为建厂地点，A1生产规模为192860吨/年，A3生产规模为96430吨/年。

对于主要原材料ST21运输，A1市的310000吨全部用于本地生产，A2市的230000吨全部运输到A1生产，A3市的270000吨也全部用于本地生产。

对于超轻型合金的销售，A1市本地销售120000吨/年，A1市运输到A2
市销售72860吨/
年，A3市生产的96430吨/年全部运输到A2销售；则A1销售量为120000吨/年，A2销售量为169290吨/年。

这样可以达到最大年平均利润2.2184*10^9元。

图3：
五、模型的评价与推广
5.1优点：
1、本文在正确、清楚地分析了题意的基础上，建立了合理、科学的线性规
划模型，便于计算得到最大年平均利润。

2、在假设基础上建立的可供计算的模型，巧妙解决了各种费用随时间的变
化等不确定问题。

4、运用了合理的数据处理方法，很好的解决了数值近似问题。

5、对于各项数据以表格或者示意图处理，给人以直观清晰的感觉，使得论
文更具说服力。

6、建立的线性规划模型与实际相联系，模型相对简单，利于操作。

结合实
际情况对问题进行求解，使得模型具有很好的通用性和推广型。

5.2缺点：
1、在假设中，我们做出了售价以及人工费等不变的假设，这与实际情况不
完全相符。

2、模型建立过程中引入的变量过多，容易引起“维数灾”，且不利于编程
处理。

3、这个模型的建立过程，许多变量受到约束不足，导致结果相对特殊，不
能适用于一般模型。

4、线性规划模型的约束条件有点简单。

六、参考文献
[1] 陈国华韦程东蒋建初付军，数学模型与数学建模方法，南开大学出版社，2012.
[2] 杨高波，MATLAB 7.0混合编程，电子工业出版社， 2006.
[3] 于润伟，MATLAB 基础及应用， 2007.
七、附录
附录1
A1,A2
distance(29+46/60+22.62/3600,53+9/60+35.81/3600,30+2/60+51.43/3600,51 +49/60+13.88/3600)*6371*pi/180
ans =
132.6625
A1,A3
distance(29+46/60+22.62/3600,53+9/60+35.81/3600,32+32/60+41.42/3600,5 2+51/60+49.18/3600)*6371*pi/180
ans =
309.5062
A2,A3
distance(30+2/60+51.43/3600,51+49/60+13.88/3600,32+32/60+41.42/3600,5 2+51/60+49.18/3600)*6371*pi/180
ans =
294.8311
附录2
tic %time begin
clear
theta=[(29+46/60+22.62/3600),(30+2/60+51.43/3600),(32+32/60+41.42/3 600),(29+46/60+22.62/3600)]/180*pi;
fai=([(53+9/60+35.81/3600),(51+49/60+13.88/3600),(52+51/60+49.18/360 0),(53+9/60+35.81/3600)]+180)/180*pi;
R=6371;
x=R*cos(theta).*cos(fai);
y=R*cos(theta).*sin(fai);
z=R*sin(theta);
for i=1:4
cor{i}=[x(i),y(i),z(i)];
end
for i=1:3
th{i}=acos(dot(cor{i},cor{i+1})/(norm(cor{i})*norm(cor{i+1})));
end
A=R*[th{1},th{2}];
distance_1={'A1到A2长度',fix(sum(A));'A1到A3长度',fix(th{3}*R)} ; %长度计算
distance_2={'A1到A2','A2到A3',fix(A(1)),fix(A(2))};
theta2=(-90:5:90)*pi/180;
fai2=(-180:5:180)*pi/180;
X=R*cos(theta2)'*cos(fai2);
Y=R*cos(theta2)'*sin(fai2);
Z=R*sin(theta2)'*ones(size(fai2));
load('topo.mat','topo','topomap1');
whos topo topomap1；
contour(0:359,-89:90,topo,[0 0],'b');
axis equal
box on
set(gca,'XLim',[0 360],'YLim',[-90 90], ...
'XTick',[0 60 120 180 240 300 360], ...
'Ytick',[-90 -60 -30 0 30 60 90]);
hold on
image([0 360],[-90 90],topo,'CDataMapping', 'scaled'); colormap(topomap1);
cla reset
axis square off
props.AmbientStrength = 0.1;
props.DiffuseStrength = 1;
props.SpecularColorReflectance =0.5;
props.SpecularExponent = 20;
props.SpecularStrength = 1;
props.FaceColor= 'texture';
props.EdgeColor = 'none';
props.FaceLighting = 'phong';
props.Cdata = topo;
surface(X,Y,Z,props); %画地球表面并进行光照和材质处理axis off
axis equal
light('position',[-8000 300 9200]);
light('position',[8000 9000 -9200]);
light('position',[800 -500 200]);
light('position',[8500 5000 -8600], 'color', [.6 .2 .2]);
view(130,150)
hold on
theta3=31/180*pi; %CD位置
fai3=103/180*pi+pi;
x3=R*cos(theta3).*cos(fai3);
y3=R*cos(theta3).*sin(fai3);
z3=R*sin(theta3);
theta4=48.194/180*pi; %LD位置
fai4=2/180*pi+pi;
x4=R*cos(theta4).*cos(fai4);
y4=R*cos(theta4).*sin(fai4);
z4=R*sin(theta4);
scatter3([x(1:3),x3,x4,0,0],[y(1:3),y3,y4,0,0],[z(1:3),z3,z4,R,-R],'filled') %画各点分布
t=(0:20)/20;
for i=1:2 %画A1到A2到A3
X=(1-t)*x(i)+t*x(i+1);
Y=(1-t)*y(i)+t*y(i+1);
Z=(1-t)*z(i)+t*z(i+1);
r=sqrt(X.*X+Y.*Y+Z.*Z);
X=R*X./r;Y=R*Y./r;Z=R*Z./r;
hold on
plot3(X,Y,Z,'r','LineWidth',2)
end
X=(1-t)*x(3)+t*x(3+1); %画A1到A3曲线
Y=(1-t)*y(3)+t*y(3+1);
Z=(1-t)*z(3)+t*z(3+1);
r=sqrt(X.*X+Y.*Y+Z.*Z);
X=R*X./r;Y=R*Y./r;Z=R*Z./r;
hold on
plot3(X,Y,Z,'b','LineWidth',2)
text(x(1),y(1),z(1),'A1') %各位置标注
text(x(2),y(2),z(2),'A2')
text(x(3),y(3),z(3),'A3')
text(x3,y3,z3,'CD')
text(0,0,R,'北极')
text(0,0,-R,'南极')
text(x4,y4,z4,'LD')
str1=num2str(fix(sum(A)));str2=num2str(fix(th{3}*R));
title(strcat('A1到A3',str1,'4公里，','A1到A3距离 ',str2,'公里'))
yu1=num2str(fix(A(1)));yu2=num2str(fix(A(2)));
text(10000,100,15000,'分段路程表示') %分段路程表示
text(9000,100,11000,strcat('A1到A2',yu1,'公里'))
text(9000,100,7000,strcat('A2到A3距离',yu2,'公里'))
format rat
distance_1
distance_2
toc
附录3
close all;
clc;
plot([(53+9/60+35.81/3600) (51+49/60+13.88/3600) (52+51/60+49.18/3600) (53+9/60+35.81/3600)],[(29+46/60+22.62/3600) (30+2/60+51.43/3600)
(32+32/60+41.42/3600) (29+46/60+22.62/3600)]);
text((53+9/60+35.81/3600),(29+46/60+22.62/3600),'A1');
text((51+49/60+13.88/3600),(30+2/60+51.43/3600),'A2');
text((52+51/60+49.18/3600),(32+32/60+41.42/3600),'A3');
axis equal
grid on
附录4
clear;
close all;
clc;
f=[-8763.126 -4781.61 -9052.6 0 636.78 1485.6297 636.78 0 1415.1892 1485.6297 1415.1892 0 0 689.845 689.845 0 1609.4322 1533.1217];
A=[2.8 0 0 0 1 1 -1 0 0 -1 0 0 0 0 0 0 0 0;0 2.8 0 0 -1 0 1 0 1 0 -1 0 0 0 0 0 0 0;0 0 2.8 0 0 -1 0 0 -1 1 1 0 0 0 0 0 0 0;1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;0 1 0 0 0 0 0 0 0 0 0 0 0 1 -1 0 0 1;1 0 0 0 0 0 0 0 0 0 0 0 0 -1 1 0 1 0];
b=[310000;230000;270000;810000/2.8;60000;170000;120000];
aeq=[0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 -1;1 0 0 0 0 0 0 0 0 0 0 0 -1 -1 0 0 0 0;0 1 0 0 0 0 0 0 0 0 0 0 0 0 -1 -1 0 0;2.8 0 0 -1 0 0 -1 0 0 -1 0 0 0 0 0 0 0 0;0 2.8 0 0 -1 0 0 -1 0 0 -1 0 0 0 0 0 0 0;0 0 2.8 0 0 -1 0 0 -1 0 0 -1 0 0 0 0 0 0];
beq=[0;0;0;0;0;0];
lb=[0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0];
[x,fval] = linprog(f,A,b,aeq,beq,lb,[])。