最优化实验的两个案例.doc

最优化实验的两个案例

1       实验案例 1.1           高速公路问题(简化)
A城和B城之间准备建一条高速公路，B城位于A城正南20公里和正东30公里交汇处，它们之间有东西走向连绵起伏的山脉。公路造价与地形特点有关，图4.2.4给出了整个地区的大致地貌情况，显示可分为三条沿东西方向的地形带。
你的任务是建立一个数学模型，在给定三种地形上每公里的建造费用的情况下，确定最便宜的路线。图中直线AB显然是路径最短的，但不一定最便宜。而路径ARSB过山地的路段最短，但是否是最好的路径呢？

平原

R · P 高地

高山

高地

S 平原

A

                                           B
图8.2  高速公路修建地段
1.1.1        问题分析
在建设高速公路时，总是希望建造费用最小。如果要建造的起点、终点在同一地貌中，那么最佳路线则是两点间连接的线段，这样费用则最省。因此本问题是一个典型的最优化问题，以建造费用最小为目标，需要做出的决策则是确定在各个地貌交界处的汇合点。
1.1.2        变量说明
：在第i个汇合点上的横坐标（以左下角为直角坐标原点），i＝1，2，…，4；x5＝30（指目的地B点的横坐标）
x=[x1，x2，x3，x4]T

li ：第i段南北方向的长度（i＝1，2，…，5）
Si：在第i段上地所建公路的长度（i＝1，2，…，5）
由问题分析可知，

C1：平原每公里的造价（单位：万元/公里）
C2：高地每公里的造价（单位：万元/公里）
C3：高山每公里的造价（单位：万元/公里）
1.1.3        模型假设
1、  假设在相同地貌中修建高速公路，建造费用与公路长度成正比；
2、  假设在相同地貌中修建高速公路在一条直线上。在理论上，可以使得建造费用最少，当然实际中一般达不到。
1.1.4        模型建立
在A城与B城之间建造一条高速公路的问题可以转化为下面的非线性规划模型。优化目标是在A城与B城之间建造高速公路的费用。

1.1.5        模型求解
这里采用Matlab编程求解。
模型求解时，分别取Ci(i=1,2,3)如下。
平原每公里的造价C1＝400万元/公里；
高地每公里的造价C2＝800万元/公里；
高山每公里的造价C3＝1200万元/公里。
输入主程序model_p97.m，运行结果如下：
model_p97
optans =
  2.2584e+004
len =
   38.9350
ans =
   12.1731   14.3323   15.6677   17.8269
求解程序见附录。
注：实际建模时必须查找资料来确定参数或者题目给定有数据）
6.模型结果及分析
通过求解可知，为了使得建造费用最小。建造地点的选择宜采取下列结果。
x1=12.1731，x2=14.3233，x3=15.6677，x4=17.8269
建造总费用为2.2584亿元。
总长度为38.9350公里。

1.1.6        求解模型的程序
（1）求解主程序
model_p97
function  x=model_p97        %数学建模教材  P97 高速公路
clear all
global C L
C=[400  800  1200];
L=[4 4 4 4 4];
x=fmincon('objfun_97',[1,1,1,1],[],[],[],[],zeros(1,4),ones(1,4)*30,'mycon_p97');
optans=objfun_97(x)
C=ones(3,1);
len = objfun_97(x)

（2）模型中描述目标函数的Matlab程序objfun_97.m
function obj=objfun_97(x)
global C L
obj=C(1)*sqrt(L(1)^2+x(1)^2) + C(2)*sqrt(L(2)^2+(x(2)-x(1))^2) + ...
   C(3)*sqrt(L(3)^2+(x(3)-x(2))^2) + ...
C(2)*sqrt(L(4)^2+(x(4)-x(3))^2)+C(1)*sqrt(L(5)^2+(...30-x(4))^2);
（3）模型中描述约束条件的Matlab函数mycon_p97.m
function [c,ceq]=mycon_p97(x)
c(1)=x(1)-x(2);
c(2)=x(2)-x(3);
c(3)=x(3)-x(4);
c(4)=x(4)-30;
ceq=[];

1       实验案例

问题侧重于线性规划和非线性规划方面的优化问题。
从这里的建模实例可以建立数学模型是最为关键和困难的一步，当看到这里建立起来的模型后，你会顿然觉得问题变得如此简单。因此，从这些实例中希望大家能够掌握建模方法，也不妨模仿这里的方法以应用到实际建模中去。

　 1.1      案例：应急设施的优化选址
问题（AMCM-86B ）里奥兰翘镇迄今还没有自己的应急设施。1986年该镇得到了建立两个应急设施的拨款，每个设施都把救护战、消防队和警察所合在一起。图（1）指出了1985年每个长方形街区发生应急事件的次数。在北边的L形状的区域是一个障碍，而在南边的长方形区域是一个右浅水池塘的公园。应急车辆驶过一条南北方向的街道平均要花15秒，而通过一条东西向的街道平均花20秒。你的任务是确定这两个应急设施的位置，使得总响应时间最少。

5
2
2
1
5
0
3
2
4
2
2
3
3
3
3
4
1

3
0
4
4
3
3

3
4
0
0
0
N

1

2
0
4
3
2
2
0
1
3
3
2
5
3
2
1
0
3
3
图（1） 1985年里奥兰翘镇每个长方街区应急事件的数目

1.1.1    问题分析
应急设施的位置
应急发生位置
应急车运行情况
总相应时间
应急设施到底修在每个街区的街角处还是可以在街道的任何地方？
对这个可以进行假设，适当简化处理。先解决应急设施在街角处的情形。对于应急设施可以在街道任何地方时，实际上可以证明：应急设施应设在街角处，才能使总响应时间最少。先在一定的假设条件下，简化问题，先解决简单情形，再处理复杂情形。
下面主要介绍穷举法在求解这个问题的应用。
1.1.2    问题假设
（1）   假设需求集中在每个街区的中心
（2）   假设应急设施位于街角处
（3）   图中给出的1985年应急次数有典型性，能够反映该街区应急事件出现的概率的大小；
（4）   应急车辆的响应时间只考虑在街道上行驶时间，其他因素（如转弯）可以忽略不计。
（5）   两个应急设施的功能完全相同。当应急事件出现时，只要从离事件发生地点最近的应急设施派出应急车辆即可。
（6）   执行任何一次应急任务的车辆都从某一个应急设施出发，完成任务后回到原设施。不出现从一个应急事件点直接到另一个应急事件发生点的情况（这是由于每个地点发生事件的概率都很小，因此两个地点同时发生事故的概率就更小，因此可以忽略这种情况。）

1.1.3    模型建立与求解
根据假设2，每个应急设施选在街角处，可能的位置至多有6×11=66个，则两个应急设施的组合就至多有C266=66×65/2=2145个，因此可以考虑用穷举法进行求解，即一一计算出每个选址方案的总响应时间，然后从中选出具有最小总响应时间的方案。
      建立直角坐标系，以该镇西北角街角为原点（即图的左下角），从北到南方向为X轴正向，从西到东为Y轴正向，并分别以南北、东西方向上的一个街区长度作为单位长，则街角的坐标（x,y）满足如下条件：0≤x≤10,0≤y≤5,x,y均为整数。
，均为整数。
而每个街区中心的坐标可定义为，其中满足如下条件：
，均为整数。
如何计算应急设施在点处到以点为街区中心的行驶时间？当不考虑障碍合池塘的影响时，行驶时间为：
如何计算应急设施在点(x,y)处到以(i+0.5,j+0.5)点为街区中心的行驶时间t(x,y,i,j)？当不考虑障碍合池塘的影响时，行驶时间为：
t(x,y,i,j)=15(|x-i-0.5|-0.5) + 20(|y-j-0.5|-0.5)
=15|x-i-0.5| + 20|y-j-0.5|-17.5

单位：秒
记为以为中心的街区的事故发生频率。如果应急设施设在两点，则总响应时间为：
记p(i,j)为以(i+0.5,j+0.5)为中心的街区的事故发生频率。如果应急设施设在(x1,y1),(x2,y2)两点，则总响应时间为：

以上模型求解不难用计算机编程实现。运用MATLAB编程得到如下结果为最优：两个应急设施应设在点（2，2），（6，2）处。

以上结论是在没有考虑障碍区域和池塘的影响的前提下得到的最优解，实际上即使考虑这两个影响，从这两个点到任何街区都可以避开L形障碍区域和池塘，所求的这两点就是最优选址。

1.1.4    更进一步结果分析
l       选址地点多余2个呢？

1.1.5    模型求解的Matlab程序
以上模型求解的Matlab程序如下：
function mymain
%应急设施的优化选址程序

p=[...
      5 2 2 1 5 0 3 2 4 2;
   2 3 3 3 3 4 1 3 0 4;
   4 3 3 0 3 4 0 0 0 0;
   1 2 0 0 4 3 2 2 0 1;
   3 3 2 5 3 2 1 0 3 3]';
%计算频数
p=p/sum(sum(p))

optvalue=inf;%初始化为无穷大
for x1=0:10,
   for y1=0:5,
      for x2=0:10,
         for y2=0:5,
            %计算每个方案的总响应时间
            curvalue=0;%当前方案的总响应时间初始化
            for i=0:9,
               for j=0:4,
                  curvalue=curvalue p(i+1,j+1)* ...
min(caltime(x1,y1,i,j),caltime(x2,y2,i,j));
               end
            end  
            if curvalue