使用MATLAB解算炼油厂的选址

背景

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

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

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

试题参考

如下图所示:

问题分析

问题一:

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

问题二:

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

问题三:

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

模型假设

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

模型建立与求解

第一问模型:

        设选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;

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

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)

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);

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 =

     5

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

            %再次生成随机数

            x(1)=82*rand;

           % x(2)=(1-0.9397)*rand+0.9397;

    end;

    while x(2)>82 | x(2)<0

        x(2)=82*rand;

    end;

    while x(3)>82 | x(3)<0

        x(3)=82*rand;

    end;

    while x(4)>82 | x(4)<0

        x(4)=82*rand;

    end;

   

    temp=myobjfunc(x);

    if temp<MIN

        MIN=temp;

        y=x;

    end;

    LIMIT=LIMIT-1;

end;

MIN

y

第三问结果:

利用蒙特卡洛算法matlab编程可求

程序多次运行的结果如下:

MIN =

  6.5582e+003

y =

   43.2099   25.0713    3.8432   81.0215

MIN =

  6.5665e+003

y =

    3.9281   81.0393   42.3370   27.9608

MIN =

  6.5583e+003

y =

    4.1848   81.1427   41.9258   25.2987

MIN =

  6.5597e+003

y =

   40.9820   25.7891    3.8891   80.8089

MIN =

  6.5713e+003

y =

   39.3252   23.5597    4.2433   80.8989

MIN =

  6.5511e+003

y =

    4.0086   81.1070   41.8500   25.6334

MIN =

  6.5553e+003

y =

   40.9798   26.6722    4.0593   80.9958

MIN =

  6.5536e+003

y =

   39.8126   24.6946    4.0317   81.0556

由多次运行程序的结果可知:

取其中最小的一次结果即可认为是本题求解答案:

MIN =6.5511e+003

y =4.0086   81.1070   41.8500   25.6334

即总费用最低为6551.1

两个炼油厂的坐标地址分别为(4.0086,81.1070),(41.8500,25.6334

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.mzph.cn/news/58152.shtml

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

curl请求https|http网站时出现Binary output can mess up your terminal

请求网站时出现​ 那么这里有几种情况 文件本身为二进制文件内容压缩 如果是第一种情况&#xff0c;那么直接保存你要下载的二进制文件&#xff0c;使用 curl https://a.com -o 文件名保存在一个文件中 或者使用 -o -直接输出在终端 curl https://a.com -o -如果你本来访问…

UE4/5的Custom节点:在VScode使用HLSL(新手入门用)

目录 custom节点 VSCode环境安装 将VSCode里面的代码放入Custom中 custom节点 可以看到这是一个简单的Custom节点&#xff1a; 而里面是可以填写代码的&#xff1a; 但是在这里面去写代码会发现十分的繁琐【按下enter后&#xff0c;不会换行&#xff0c;也不会自动缩进】 …

火山引擎发布自研视频编解码芯片

2023年8月22日&#xff0c;火山引擎视频云宣布其自研的视频编解码芯片已成功出片。经验证&#xff0c;该芯片的视频压缩效率相比行业主流硬件编码器可提升30%以上&#xff0c;未来将服务于抖音、西瓜视频等视频业务&#xff0c;并将通过火山引擎视频云开放给企业客户。 火山引…

【网络】多路转接——五种IO模型 | select

&#x1f431;作者&#xff1a;一只大喵咪1201 &#x1f431;专栏&#xff1a;《网络》 &#x1f525;格言&#xff1a;你只管努力&#xff0c;剩下的交给时间&#xff01; 五种IO模型 | select &#x1f367;五种IO模型&#x1f367;select&#x1f9c1;认识接口&#x1f9c1…

视频中的声音怎么提取出来?这样做提取出来很简单

提取视频中的声音可以有多种用途。例如&#xff0c;我们可能希望从视频中提取音乐或音效&#xff0c;以在其他项目中使用。或者&#xff0c;可能需要将视频中的对话转录为文本&#xff0c;以便更轻松地编辑和共享内容。无论目的是什么&#xff0c;提取视频中的声音都可以帮助我…

调用自实现MyGetProcAddress获得CreateFileA函数并调用创建写入文件

写文件如下 #include <iostream> #include <Windows.h>typedef HANDLE(WINAPI* CreateFileAFunc)(LPCSTR, DWORD, DWORD, LPSECURITY_ATTRIBUTES, DWORD, DWORD, HANDLE);DWORD MyGetProcAddress(_In_ HMODULE hModule,_In_ LPCSTR lpProcName ){PIMAGE_DOS_HEADE…

Mycat教程+面试+linux搭建

目录 一 MyCAT介绍 二 常见的面试题总结 三 linux下搭建Mycat 一 MyCAT介绍 1.1. 什么是MyCAT&#xff1f; 简单的说&#xff0c;MyCAT就是&#xff1a; 一个彻底开源的&#xff0c;面向企业应用开发的“大数据库集群” 支持事务、ACID、可以替代Mysql的加强版数据库 一个可…

uni-app里使用webscoket

实现思路和vue中是一样的。如果想看思路可以看这篇文章&#xff1a;websocket 直接上可以运行的代码&#xff1a; 一、后端nodeJS代码&#xff1a; 1、新建项目文件夹 2、初始化项目&#xff1a; npm init -y 3、项目里安装ws npm i ws --save 4、nodeJS代码&#xff1…

SmartInspect Professional .Net Delphi Crack

SmartInspect Professional .Net & Delphi Crack SmartInspect Professional是一个用于调试和跟踪.NET、Java和Delphi软件的高级日志记录工具。它使您能够识别错误&#xff0c;找到客户问题的解决方案&#xff0c;并让您清楚地了解软件在不同环境和条件下的工作方式。可以轻…

Redis 7 第三讲 数据类型 进阶篇

⑥ *位图 bitmap 1. 理论 由0和1 状态表现的二进制位的bit 数组。 说明:用String 类型作为底层数据结构实现的一种统计二值状态的数据类型 位图本质是数组,它是基于String 数据类型的按位操作。该数组由多个二进制位组成,每个二进制位都对应一个偏…

3、监测数据采集物联网应用开发步骤(3)

监测数据采集物联网应用开发步骤(2) 系统整体结构搭建 新建项目 输入项目名称&#xff1a;MonitorData 所谓兵马未动粮草先行&#xff0c;按下图创建好对应的模块备用&#xff1a; com.plugins 业务插件模块 com.zxy.adminlog 日志或文本文…

基于python+pyqt的opencv汽车分割系统

目录 一、实现和完整UI视频效果展示 主界面&#xff1a; 识别结果界面&#xff1a; 查看分割处理过程图片界面&#xff1a; 二、原理介绍&#xff1a; 加权灰度化 ​编辑 二值化 滤波降噪处理 锐化处理 边缘特征提取 图像分割 完整演示视频&#xff1a; 完整代码链…

计算机竞赛 基于YOLO实现的口罩佩戴检测 - python opemcv 深度学习

文章目录 0 前言1 课题介绍2 算法原理2.1 算法简介2.2 网络架构 3 关键代码4 数据集4.1 安装4.2 打开4.3 选择yolo标注格式4.4 打标签4.5 保存 5 训练6 实现效果6.1 pyqt实现简单GUI6.3 视频识别效果6.4 摄像头实时识别 7 最后 0 前言 &#x1f525; 优质竞赛项目系列&#xf…

【推荐】Spring与Mybatis集成整合

目录 1.概述 2.集成 2.1代码演示&#xff1a; 3.整合 3.1概述 3.2 进行整合分页 接着上两篇&#xff0c;我已经写了Mybatis动态之灵活使用&#xff0c;mybatis的分页和特殊字符的使用方式接下来把它们集成起来&#xff0c;是如何的呢&#x1f447;&#x1f447;&#x1…

vscode 无法跳转第三方安装包

vscode 无法跳转第三方安装包 场景&#xff1a;使用vscode写代码时&#xff0c; 第三方的安装包无法使用ctrl 左键&#xff0c;点击进入查看&#xff0c; 不方便源码查看 解决办法&#xff1a; 使用快捷键 Ctrl Shift P&#xff0c; 进入命令搜索框搜索 setting.json 编辑…

react解决死循环方法?

使用useeffect&#xff08;副作用&#xff09;方法结束这个操作 1、导入useeffect、useState 2、把下方代码写入&#xff1a;里面填写的是你要终止某个东西的代码 注意&#xff1a;不可不写&#xff0c;也可以写依赖或不写

精进面试技巧:如何在程序员面试中脱颖而出

&#x1f337;&#x1f341; 博主猫头虎 带您 Go to New World.✨&#x1f341; &#x1f984; 博客首页——猫头虎的博客&#x1f390; &#x1f433;《面试题大全专栏》 文章图文并茂&#x1f995;生动形象&#x1f996;简单易学&#xff01;欢迎大家来踩踩~&#x1f33a; &a…

《人月神话》:chapter 4 系统设计中的“专制”和“民主”

以下总结来自于《人月神话》 第四章 &#xff1a;贵族制&#xff0c;民主制和系统设计 系统设计中最重要的因素&#xff1a;概念完整性 1.设计必须由一个人或者具有共识的小型团队来完成 2.大型系统架构设计与具体实现相分离 3.必须有人控制概念&#xff0c;确保完整性&…

Verilog基础:块语句

相关阅读 Verilog基础专栏https://blog.csdn.net/weixin_45791458/category_12263729.html?spm1001.2014.3001.5482 1、块语句 块语句(block statements)是一种把语句组织在一起&#xff0c;这样他们在语法上就像单个语句一样工作。Verilog HDL中有两种类型的块&#xff1a; …

树状表格子节点移动 - 在Vue.js中实现上下移动子节点的表格功能

目录 功能介绍 示例 代码 视图部分 逻辑部分 完整代码 功能介绍 本文介绍了如何在Vue.js框架下实现一个树状表格&#xff0c;其中支持选择子节点行的上下移动。通过这个功能&#xff0c;用户可以方便地改变子节点的顺序。代码示例和详细的实现步骤将展示如何使用Vue.js的相…