背景

记得有一年的数据建模大赛，试题是炼油厂的选址，最后我们采用MATLAB编写（复制）蒙特卡洛算法，还到了省级一等奖，这里把仅有一些记忆和材料，放到这里来，用来纪念消失的青春。

本文使用素材下载，内含MATLAB代码

使用蒙特卡洛算法解算炼油厂的选址MATLAB程序，提供试题照片，以及MATLAB代码资源-CSDN文库

试题参考

如下图所示：

问题分析

问题一：

本问的炼油厂选址是九口油井的任一处，我们可以把九口油井依次作为炼油厂，然后分别计算其费用。

问题二：

本问的炼油厂选址范围应在0<横坐标x<100、0<纵坐标y<100,的矩形区域内。可以把问题转化为：在该矩形区域内找一点使其满足总费用最少。

问题三：

本问的两个炼油厂的选址纵横坐标都应在[0,82]的范围内。问即可转化成在该矩形区域内找两点，使其总费用最低。

模型假设

总费用与距离成正比，设比例系数为1；
总费用与油量成正比，设比例系数为1；
其他因素不影响总费用；
一口油井的原油智能运往一个炼油厂；

模型建立与求解

第一问模型：

设选1号油井为炼油厂产生的费用为feiyoug1；

选1号油井为炼油厂产生的费用为feiyoug1；

选2号油井为炼油厂产生的费用为feiyoug2；

选3号油井为炼油厂产生的费用为feiyoug3；

选4号油井为炼油厂产生的费用为feiyoug4；

…

选9号油井为炼油厂产生的费用为feiyoug9；

由问题一分析可知问题可转化为求

feiyong1=(abs(x2-x1)+abs(y2-y1))*chanliang2+(abs(x3-x1)+abs(y3-y1))*chanliang3+(abs(x4-x1)+abs(y4-y1))*chanliang4+(abs(x5-x1)+abs(y5-y1))*chanliang5+(abs(x6-x1)+abs(y6-y1))*chanliang6+(abs(x7-x1)+abs(y7-y1))*chanliang7+(abs(x8-x1)+abs(y8-y1))*chanliang8+(abs(x9-x1)+abs(y9-y1))*chanliang9

（其中abs是求绝对值，其余八个油井到1号油井的折线距离与产量之积作为费用，总费用即为八个油井产生费用之和）

同法可求feiyong2、feiyong3、……、feiyong9。

第一问求解：

利用matlab编程求解如下：

第一问程序：

x1=22;

y1=38;

x2=8;

y2=13;

x3=4;

y3=81;

x4=51;

y4=32;

x5=38;

y5=11;

x6=17;

y6=12;

x7=81;

y7=63;

x8=19;

y8=45;

x9=62;

y9=12;

chanliang1=17;

chanliang2=40;

chanliang3=60;

chanliang4=20;

chanliang5=25;

chanliang6=15;

chanliang7=50;

chanliang8=8;

chanliang9=30;

feiyong2=(abs(x1-x2)+abs(y1-y2))*chanliang1+(abs(x3-x2)+abs(y3-y2))*chanliang3+(abs(x4-x2)+abs(y4-y2))*chanliang4+(abs(x5-x2)+abs(y5-y2))*chanliang5+(abs(x6-x2)+abs(y6-y2))*chanliang6+(abs(x7-x2)+abs(y7-y2))*chanliang7+(abs(x8-x2)+abs(y8-y2))*chanliang8+(abs(x9-x2)+abs(y9-y2))*chanliang9

feiyong3=(abs(x1-x3)+abs(y1-y3))*chanliang1+(abs(x2-x3)+abs(y2-y3))*chanliang2+(abs(x4-x3)+abs(y4-y3))*chanliang4+(abs(x5-x3)+abs(y5-y3))*chanliang5+(abs(x6-x3)+abs(y6-y3))*chanliang6+(abs(x7-x3)+abs(y7-y3))*chanliang7+(abs(x8-x3)+abs(y8-y3))*chanliang8+(abs(x9-x3)+abs(y9-y3))*chanliang9

feiyong4=(abs(x1-x4)+abs(y1-y4))*chanliang1+(abs(x2-x4)+abs(y2-y4))*chanliang2+(abs(x3-x4)+abs(y3-y4))*chanliang3+(abs(x5-x4)+abs(y5-y4))*chanliang5+(abs(x6-x4)+abs(y6-y4))*chanliang6+(abs(x7-x4)+abs(y7-y4))*chanliang7+(abs(x8-x4)+abs(y8-y4))*chanliang8+(abs(x9-x4)+abs(y9-y4))*chanliang9

feiyong5=(abs(x1-x5)+abs(y1-y5))*chanliang1+(abs(x2-x5)+abs(y2-y5))*chanliang2+(abs(x3-x5)+abs(y3-y5))*chanliang3+(abs(x4-x5)+abs(y4-y5))*chanliang4+(abs(x6-x5)+abs(y6-y5))*chanliang6+(abs(x7-x5)+abs(y7-y5))*chanliang7+(abs(x8-x5)+abs(y8-y5))*chanliang8+(abs(x9-x5)+abs(y9-y5))*chanliang9

feiyong6=(abs(x1-x6)+abs(y1-y6))*chanliang1+(abs(x2-x6)+abs(y2-y6))*chanliang2+(abs(x3-x6)+abs(y3-y6))*chanliang3+(abs(x4-x6)+abs(y4-y6))*chanliang4+(abs(x5-x6)+abs(y5-y6))*chanliang5+(abs(x7-x6)+abs(y7-y6))*chanliang7+(abs(x8-x6)+abs(y8-y6))*chanliang8+(abs(x9-x6)+abs(y9-y6))*chanliang9

feiyong7=(abs(x1-x7)+abs(y1-y7))*chanliang1+(abs(x2-x7)+abs(y2-y7))*chanliang2+(abs(x3-x7)+abs(y3-y7))*chanliang3+(abs(x4-x7)+abs(y4-y7))*chanliang4+(abs(x5-x7)+abs(y5-y7))*chanliang5+(abs(x6-x7)+abs(y6-y7))*chanliang6+(abs(x8-x7)+abs(y8-y7))*chanliang8+(abs(x9-x7)+abs(y9-y7))*chanliang9

feiyong8=(abs(x1-x8)+abs(y1-y8))*chanliang1+(abs(x2-x8)+abs(y2-y8))*chanliang2+(abs(x3-x8)+abs(y3-y8))*chanliang3+(abs(x4-x8)+abs(y4-y8))*chanliang4+(abs(x5-x8)+abs(y5-y8))*chanliang5+(abs(x6-x8)+abs(y6-y8))*chanliang6+(abs(x7-x8)+abs(y7-y8))*chanliang7+(abs(x9-x8)+abs(y9-y8))*chanliang9

feiyong9=(abs(x1-x9)+abs(y1-y9))*chanliang1+(abs(x2-x9)+abs(y2-y9))*chanliang2+(abs(x3-x9)+abs(y3-y9))*chanliang3+(abs(x4-x9)+abs(y4-y9))*chanliang4+(abs(x5-x9)+abs(y5-y9))*chanliang5+(abs(x6-x9)+abs(y6-y9))*chanliang6+(abs(x7-x9)+abs(y7-y9))*chanliang7+(abs(x8-x9)+abs(y8-y9))*chanliang8

第一问运行结果：

feiyong1 =

13720

feiyong2 =

15317

feiyong3 =

18635

feiyong4 =

14835

feiyong5 =

15185

feiyong6 =

14857

feiyong7 =

20108

feiyong8 =

13980

feiyong9 =

16970

由程序运行结果可知：

feiyong1最小，即炼油厂建在1号油井附近，总费用最少。

第二问模型：

由分析可设该炼油厂得地址为（x(1),x(2)）,且0<x(1)<100,0<x(2)<100

即求zongfeiyong最小

zongfeiyong=17*sqrt((x(1)-22)^2+(x(2)-38)^2)+40*sqrt((x(1)-8)^2+(x(2)-13)^2)+60*sqrt((x(1)-4)^2+(x(2)-81)^2)+20*sqrt((x(1)-51)^2+(x(2)-32)^2)+25*sqrt((x(1)-38)^2+(x(2)-11)^2)+15*sqrt((x(1)-17)^2+(x(2)-12)^2)+50*sqrt((x(1)-81)^2+(x(2)-63)^2)+8*sqrt((x(1)-19)^2+(x(2)-45)^2)+30*sqrt((x(1)-62)^2+(x(2)-12)^2);

其中sqrt为平方的意思，分别求的每个油井的到炼油厂的距离与其产量的积，再作和，即为总费用。

第二问用matlab编程求解如下（利用matlab总求最优解函数fmincon）：

第二问程序：

function zongfeiyong=timu2(x)

xmin=[0;0];

xmax=[100;100];

[x,fopt,flag,c]=fmincon('timu2',zeros(2,1),[],[],[],[],xmin,xmax)

第二问运行结果：

x =

32.4224

35.0597

fopt =

1.0213e+004

flag =

c =

iterations: 8

funcCount: 26

lssteplength: 1

stepsize: 2.0001e-004

algorithm: [1x44 char]

firstorderopt: 5.9487e-004

constrviolation: -32.4224

message: [1x844 char]

由程序运行结果可知：

即选址坐标为（32.4224，35.0597）

此时总费用最少为10213

第三问模型：

我们可以在纵、横坐标都为[0,82]的范围内任意选取两点做炼油厂。总费用即可转换为RES=min(17*sqrt((x(1)-22)^2+(x(2)-38)^2),17*sqrt((x(3)-22)^2+(x(4)-38)^2))+min(40*sqrt((x(1)-8)^2+(x(2)-13)^2),40*sqrt((x(3)-8)^2+(x(4)-13)^2))+min(60*sqrt((x(1)-4)^2+(x(2)-81)^2),60*sqrt((x(3)-4)^2+(x(4)-81)^2))+min(20*sqrt((x(1)-51)^2+(x(2)-32)^2),20*sqrt((x(3)-51)^2+(x(4)-32)^2))+min(25*sqrt((x(1)-38)^2+(x(2)-11)^2),25*sqrt((x(3)-38)^2+(x(4)-11)^2))+min(15*sqrt((x(1)-17)^2+(x(2)-12)^2),15*sqrt((x(3)-17)^2+(x(4)-12)^2))+min(50*sqrt((x(1)-81)^2+(x(2)-63)^2),50*sqrt((x(3)-81)^2+(x(4)-63)^2))+min(8*sqrt((x(1)-19)^2+(x(2)-45)^2),8*sqrt((x(3)-19)^2+(x(4)-45)^2))+min(30*sqrt((x(1)-62)^2+(x(2)-12)^2),30*sqrt((x(3)-62)^2+(x(4)-12)^2))

其中min（A,B）是求A、B中较小的数。

第三问求解：

利用蒙特卡洛算法，我们对x(1),x(2),x(3), x(4)随机生成一千万组（范围在[0,82]内），带入总费用公式RES中。求得其中最小值。利用matlab编程如下（程序多次运行）：

第三问程序：

% 原函数

function RES=myobjfunc(x)

RES=min(17*sqrt((x(1)-22)^2+(x(2)-38)^2),17*sqrt((x(3)-22)^2+(x(4)-38)^2))+min(40*sqrt((x(1)-8)^2+(x(2)-13)^2),40*sqrt((x(3)-8)^2+(x(4)-13)^2))+min(60*sqrt((x(1)-4)^2+(x(2)-81)^2),60*sqrt((x(3)-4)^2+(x(4)-81)^2))+min(20*sqrt((x(1)-51)^2+(x(2)-32)^2),20*sqrt((x(3)-51)^2+(x(4)-32)^2))+min(25*sqrt((x(1)-38)^2+(x(2)-11)^2),25*sqrt((x(3)-38)^2+(x(4)-11)^2))+min(15*sqrt((x(1)-17)^2+(x(2)-12)^2),15*sqrt((x(3)-17)^2+(x(4)-12)^2))+min(50*sqrt((x(1)-81)^2+(x(2)-63)^2),50*sqrt((x(3)-81)^2+(x(4)-63)^2))+min(8*sqrt((x(1)-19)^2+(x(2)-45)^2),8*sqrt((x(3)-19)^2+(x(4)-45)^2))+min(30*sqrt((x(1)-62)^2+(x(2)-12)^2),30*sqrt((x(3)-62)^2+(x(4)-12)^2));

%蒙特卡罗算法

MIN=inf;

%取一千万组随机数

LIMIT=10000000;

while LIMIT>0

%生成指定范围内的随机数

x(1)=82*rand;

x(2)=82*rand;

x(3)=82*rand;

x(4)=82*rand;

%过滤不符合条件的随机数

while x(1)>82 | x(1)<0 %| x(2)>81 | x(2)<0 | x