数学实验约束规划.docx

实验7约束优化
分1黄浩
一、 实验目的
1. 学会用MATLAB工具箱求解非线性规划的方法。
2. 练习建立实际问题的非线性规划模型。
二、 实验内容
1. 《数学实验》第一版（问题5）
问题叙述：
对问题：
min100(x2-x12)2+(1-x1)2+90(x4-x32)2+(1-x3)2+10.11-x22+1-x42+19.8x2-1(x4-1)，增长以下条件，并分别取初值(-3，-1，-3，-1)和(3，1，3，1)，求解非线性优化：
(1)-10≤xi≤10；
(2)-10≤xi≤10，x1x2-x1-x2+1.5≤0，x1x2+10≥0，-100≤x1x2x3x4≤100；
(3)-10≤xi≤10，x1x2-x1-x2+1.5≤0，x1x2+10≥0，x1+x2=0，x1x2x3x4=16。
再试取不同的初值或用分析梯度计算，比较计算结果，你能从中得到什么启示？
实验过程：
（1）小题
先编写带有梯度的目的函数（程序见四.1），分别设初值为(-3,-1,-3,-1)、(3,1,3,1)和(6,-3,9,-12)，使用数值方法和分析方法分别求解（程序见四.2），所得结果如下：
初值
方法
最优解(x1,x2,x3,x4)
最优值z
迭代次数
调用次数
(-3,-1,-3,-1)
数值
1.0001631
1.0003534
0.9997828
0.9995948
3.0895E-07
24
136
分析
1.0001633
1.0003538
0.9997827
0.9995947
3.08745E-07
24
63
(3,1,3,1)
数值
1.0002399
1.0005064
0.9997450
0.9995082
3.24271E-07
71
379
分析
1.0000609
1.0000821
0.9999463
0.9998689
2.42479E-07
65
153
(6,-3,9,-12)
数值
0.9999935

0.9999963
0.9999798
0.9999634
2.66218E-08
71
389
分析
1.0001596
1.0003463
0.9998086
0.9996536
3.02103E-07
73
178
由上表可见，目的函数的最优值近似为0，最优解近似为(1,1,1,1)，并且不管选取哪个初值、哪种计算方法，最终都能得到近似对的的结果。此外，由于在命令中没有对误差上限进行设立（即默认值），因此这六个结果的误差数量级都是相同的。
若进行更进一步地分析，不难发现，在相同初值和相同精度 在上一段中有解释，即精度TolFun和TolX都为相同的默认值

条件下，数值方法和分析方法的求解误差和迭代次数相差不大，但分析方法的目的函数调用次数明显少于数值方法。另一方面，在相同计算方法的条件下，不同初值的迭代次数和调用次数也不同，并且这与初值x0和最优解x之间的距离没有本质的联系，如初值为(-3,-1,-3,-1)时的迭代次数少于(3,1,3,1)，但显然前者与最优解(1,1,1,1)之间的距离要大于后者，因此，迭代次数不能单纯地用距离来衡量。
（2）小题
这一小问相对于前者增长了非线性的不等式约束，因此要增长一个非线性约束的M文献（程序见四.3），分别设初值为(-3,-1,-3,-1)、(3,1,3,1)和(1.1,1.1,1.1,1.1)，使用数值方法和分析方法分别求解（程序见四.4），所得结果如下：
初值
方法
最优解(x1,x2,x3,x4)
最优值z
迭代次数
调用次数
(-3,-1,-3,-1)
数值
1.4205956
-0.1887903
-1.4394174
2.0813389
493.799
18
111
分析
1.4205956
-0.1887903
-1.4394178
2.0813398
493.799
18
56
(3,1,3,1)
数值
1.4202330
-0.1903543
1.4660517
2.1510748
487.988
16
97
分析
1.4202330
-0.1903543
1.4660519
2.1510754
487.988
16
48
(1.1,1.1,1.1,1.1)
数值
-1.1082373
1.2371649
0.8803495
0.7742603
4.490
47
260
分析
-1.1082327
1.2371655
0.8803625
0.7742826
4.490
48
116
由于增长了新的约束条件，最优解和最优值都与上一问有很大出入，并且选取不同的初值也会导致最优解的不同，这是由于matlab只能根据一定的初值和步长求解局部最优点而导致的，获得全局最优点需要枚举所有的局部最优点，并进行比较才干得到。与上一问相似的地方是，数值方法和分析方法在解的精度和迭代次数上相近，但分析方法的目的函数调用次数较少。
（3）小题
这一小问相对于前者增长了非线性的等式约束和线性的等式约束，因此要修改约束函数的M文献（程序见四.5），分别设初值为(-3,-1,-3,-1)、(3,1,3,1)和(1.1,1.1,1.1,1.1)，使用数值方法和分析方法分别求解（程序见四.6），所得结果如下：
初值
方法
最优解(x1,x2,x3,x4)
最优值z
迭代次数
调用次数
(-3,-1,-3,-1)
数值
1.7747547
-1.7747547
1.4276319
-3.5581729
5782.576
13
77
分析
1.7747547
-1.7747547
1.4276319
-3.5581729
5782.576
13
37
(3,1,3,1)
数值
-2.2377236
2.2377236
1.1946170
-2.6747205
2353.689
12
66
分析
-2.2377237
2.2377237
1.1946170
-2.6747205
2353.689
12
29
(1.1,1.1,1.1,1.1)
数值
0.1730997

0.1730997
17.1669721
17.1669721
6935104
1075
20231
分析
0.2674696
0.2683346
10.5872573
10.5872573
928144.11
1332
20231
由上表可知，更换约束条件后，题目给定的两个初值(-3,-1,-3,-1)和(3,1,3,1)产生了不同的最优解和最优值，但都能以较少的迭代次数收敛到各自的局部最优解，并且数值方法和分析方法的精度十分接近，只是后者的目的函数调用次数较少。
然而，若使用初值(1.1,1.1,1.1,1.1)，则在目的函数调用次数上限MaxFun=20230的条件下无法完毕求解，并且搜索过程是发散的，无法达成最优解。
得出结论：
通过对这个非线性规划问题的求解，我获得了以下几点启示：
1) Matlab数值方法的求解精度很高，一般情况下没有必要使用分析方法求解梯度，后者只是在目的函数的调用次数上较少，对于精度的奉献比较小，在中小规模运算和迭代次数较少时可直接使用数值方法。
2) 在最优解唯一的情况下，使用不同初值（规定收敛）所得最优解和最优值近似，但迭代次数有较大差异，并且迭代次数与初值x0-最优解x之间的距离没有本质联系。而在最优解不唯一的情况下，不同初值也许会产生不同的最优解，它们都是局部极值点，需要比较所有局部极值点后才干得到全局极值点（最值点）。
3) 在优化实验中，对优化函数的控制参数进行设立看似多余，事实上很有必要，例如设立最大调用次数和最大迭代次数来防止发散情况下的无限迭代，以及调节输出精度来满足更高的规定等等。
4) 这次进行的是小规模的非线性优化实验，而在更深层次的工程领域，也许会出现大规模的实验，求解速度会十分缓慢，这时，使用更好的算法和更合适的初值就变得尤为重要。
2.《数学实验》第一版（例题：供应与选址）
问题描述：
某公司有6个建筑工地要开工，每个工地的位置(用平面坐标x，y表达，距离单位：km)及水泥日用量d(t)由下表给出。目前有两个临时料场位于A5，1，B(2，7)，日储量各有20t，假设从料场到工地之间均有直线道路相连，试制订天天的供应计划，即从两料场分别向各工地运送多少吨水泥，使总的吨公里数最小。
为进一步减少吨公里数，打算舍弃两个临时料场，改建两个新的，日储量仍各为20吨，问应建在何处？节省的吨公里数有多大？
工地的位置(x，y)及水泥日用量d
工地
1
2
3
4
5
6
x
1.25
8.75
0.5
5.75
3
7.25
y
1.25
0.75
4.75
5
6.5
7.75
d
3
5
4
7
6
11
模型转换及实验过程：
设：各工地的位置分别为(ai,bi)，水泥的日用量分别为di（i=1,…,6，分别表达六个工地），料场位置设为(xj,yj)，日储量分别为ej（j=1,2，分别表达A,B料场），从料场j运往工地i的水泥运送量为cij，则原问题可以表达为以下的优化问题：
minf=j=12i=16cijxj-ai2+yj-bi2
s.t.i=16cij≤ej,j=1,2
j=12cij=di,i=1,2,…,6
cij≥0,i=1,2,……,6,j=1,2
以cij作为决策变量，不难发现这是一个线性规划问题，使用matlab求解这个线性规划（程序见四.7），结果如下：
工地i
1
2
3
4
5
6
料场A
2.
4.
0.
6.
0.

0.
料场B
0.
0.
3.
0.
5.
10.
四舍五入后得到：
工地i
1
2
3
4
5
6
最优值z
料场A
3
5
0
7
0
1
136.228
料场B
0
0
4
0
6
10
因此，当各料场运送量按上表执行时，总的吨公里数为136.228，且为全局最小值。
假如舍弃现有的A、B料场，改建新料场，则现在的优化问题变成了非线性优化问题，目的函数变成了非线性的，决策变量是料场向工地的运送吨数和料场位置，而约束条件及上下界都不变。
先编写目的函数的M文献（程序见四.8），将原A、B料场的位置作为新优化问题的初值，然后使用matlab求解这个非线性约束（程序见四.9），结果如下（通过了四舍五入）：
工地i
1
2
3
4
5
6
料场位置
最优值z
料场A
3
5
4
7
1
0
(5.696,4.929)
88.883
料场B
0
0
0
0
5
11
(7.250,7.750)
在这种情况的下的总吨公里数为88.883，迭代次数为120次。
由于matlab的非线性优化只能通过一个初值来求得局部极小值，因此我们换用其他的位置初值进行求解（运送量初值均设为0,程序略），结果如下：
位置初值
工地i
1
2
3
4
5
6
料场位置
最优值z
(6,5),(4,3)
料场A
0
5
0
0
0
11
(7.250,7.750)
85.266
料场B
3
0
4
7
6
0
(3.255,5.652)
(4,4),(4,4)
料场A
1.5
2.5
2
3.5
3
5.5
(4,4)
139.613
料场B
2.5
1.5
2
3.5
3
5.5
(4,4)
(1,4),(6,9)
料场A
3
0
4
3
6
0
(2.353,5.368)
87.001
料场B
0
5
0
4
0
11
(7.250,7.750)
由上表可见，选取不同的初值都可以收敛到一定的局部最优解，然而它们的目的函数最优值却各不相同，在当前所有的局部最优解中，最小的为85.266，即上表中的第一个初值所得结果。
与此同时，我们发现j=12cij=di,i=1,2,…,6是六个线性的等式约束，因此一定可以通过消元法消除六个变量，恰好这六个等式约束都为二元变量，消元就更加容易了。重新对这个非线性规划问题进行整理，得：
mini=16{cix1-ai2+y1-bi2+(di-ci)x2-ai2+y2-bi2}
s.t.0≤ci≤di,i=1,2,……,6
i=16ci≤e1
i=16di-ci≤e2
因此，整理后的决策变量为新料场A’向六个工地的运送量ci和两个新料场的位置(xj,yj)，重新编译目的函数的M文献（程序见四.10），使用原料场的位置为位置初值，然后用matlab求解这个非线性规划（程序见四.11），结果如下：
工地i
1
2
3
4
5
6
料场位置
最优值z
料场A
3
0
4
7
6
0
(3.255,5,652)
85.266
料场B
0
5
0
0
0
11
(7.250,7.750)
由上表可见，在重新整理了这个优化模型之后，得到了更优化的解，并且发现这种情况下的迭代次数为36次，远小于一开始的120次，此外，通过与前面获得的局部最优值比较可得，85.266是当前获得的全局最优值。这说明，减少问题的维数对于提高解精度以及减小迭代次数方面有好处的。当然，这也要视情况而定，当等式约束皆为非线性约束、或是多元的线性约束时，盲目进行消元降维会增长人的计算复杂度，编程时也容易犯错，因此，我认为，只有像本题这种线性的二元约束，才有消元的必要，其他情况下盲目消元是无益的，毕竟matlab的数值功能还是很强大的。
得出结论：
通过实验可得，若使用临时料场A、B，则该问题为线性规划问题，当料场向各个工地的运送量如下表所示时，总的吨公里数最少：
工地i
1
2
3
4
5
6
最优值z
料场A
3
5
0
7

0
1
136.228
料场B
0
0
4
0
6
10
若改建新料场，则当料场位置和运送量如下表所示时，总的吨公里数最少：
工地i
1
2
3
4
5
6
料场位置
最优值z
料场A
3
0
4
7
6
0
(3.255,5,652)
85.266
料场B
0
5
0
0
0
11
(7.250,7.750)
并且通过这两个表的比较可见，新料场的吨公里数比临时料场的吨公里数减少了50.962。
三、 实验总结
本次实验是使用matlab求解带约束的非线性优化问题。具体的实行方法和之前的无约束优化、线性优化十分相似从matlab的编程角度上看，本次实验较为简朴，只需要注意程序的规范性、注意角标的对的使用即可。
当然，通过本次实验过程，我也从非线性规划的算法中获得了一些启示：一是要合理选取初值，以减少迭代次数并获得更可靠的局部最优值；二是在小规模和低迭代的条件下，使用分析方法求解梯度对于提高精度无益；三是要充足运用matlab提供的控制参数，设立精度下限和迭代上限，并对算法做出一个合理的选择；四是在约束条件中，若有二元的等式约束，则可通过消元降维来提高解的精度，并可减少迭代时间。
四、 程序清单
1. 第一题——（1）小题的目的函数
function[f,g]=pro71(x)
f=100*(x(2)-x(1)^2)^2+(1-x(1))^2+90*(x(4)-x(3)^2)^2+(1-x(3))^2+10.1*((1-x(2))^2+(1-x(4))^2)+19.8*(x(2)-1)*(x(4)-1);
ifnargout>1
g(1)=-400*x(1)*(x(2)-x(1)^2)+2*(x(1)-1);
g(2)=200*(x(2)-x(1)^2)+20.2*(x(2)-1)+19.8*(x(4)-1);
g(3)=-360*x(3)*(x(4)-x(3)^2)+2*(x(3)-1);
g(4)=180*(x(4)-x(3)^2)+20.2*(x(4)-1)+19.8*(x(2)-1);
end
end
2. 第一题——（1）小题的求解
v1=[-10,-10,-10,-10];
v2=[10,10,10,10];
opt1=optimset('largescale','off','MaxIter',3000,'MaxFun',20230);
opt2=optimset(opt1,'GradObj','on','GradConstr','on','DerivativeCheck','on');
x0=[-3,-1,-3,-1];
[x1,fv1,ef,out1]=fmincon(@pro71,x0,[],[],[],[],v1,v2,[],opt1)
x0=[-3,-1,-3,-1];
[x2,fv2,ef,out2]=fmincon(@pro71,x0,[],[],[],[],v1,v2,[],opt2)
x0=[3,1,3,1];
[x3,fv3,ef,out3]=fmincon(@pro71,x0,[],[],[],[],v1,v2,[],opt1)
x0=[3,1,3,1];
[x4,fv4,ef,out4]=fmincon(@pro71,x0,[],[],[],[],v1,v2,[],opt2)
x0=[6,-3,9,-12];
[x5,fv5,ef,out5]=fmincon(@pro71,x0,[],[],[],[],v1,v2,[],opt1)
x0=[6,-3,9,-12];
[x6,fv6,ef,out6]=fmincon(@pro71,x0,[],[],[],[],v1,v2,[],opt2)
[x1;x2;x3;x4;x5;x6]
[fv1;fv2;fv3;fv4;fv5;fv6]
[out1.iterations;out2.iterations;out3.iterations;out4.iterations;out5.iterations;out6.iterations]
[out1.funcCount;out2.funcCount;out3.funcCount;out4.funcCount;out5.funcCount;out6.funcCount]
3. 第一题——（2）小题的非线性约束函数
function[c,ceq,g,geq]=pro71con(x)
c=[x(1)*x(2)-x(1)-x(2)+1.5;-x(1)*x(2)-10;x(1)*x(2)*x(3)*x(4)-100;-x(1)*x(2)*x(3)*x(4)-100];
ceq=[0,0,0,0];
ifnargout>2
g=[x(2)-1,-x(2),x(2)*x(3)*x(4),-x(2)*x(3)*x(4);x(1)-1,-x(1),x(1)*x(3)*x4),-x(1)*x(3)*x(4);0,0,x(1)*x(2)*x(4),-x(1)*x(2)*x(4);0,0,x(1)*x(2)*x(3),-x(1)*x(2)*x(3)];
geq=zeros(4,4);
end
end
4. 第一题——（2）小题的求解
v1=[-10,-10,-10,-10];
v2=[10,10,10,10];
opt1=optimset('largescale','off','MaxIter',3000,'MaxFun',20230);
opt2=optimset(opt1,'GradObj','on','GradConstr','on','DerivativeCheck','on');
x0=[-3,-1,-3,-1];
[x1,fv1,ef,out1]=fmincon(@pro71,x0,[],[],[],[],v1,v2,@pro71con,opt1)
x0=[-3,-1,-3,-1];
[x2,fv2,ef,out2]=fmincon(@pro71,x0,[],[],[],[],v1,v2,@pro71con,opt2)
x0=[3,1,3,1];
[x3,fv3,ef,out3]=fmincon(@pro71,x0,[],[],[],[],v1,v2,@pro71con,opt1)
x0=[3,1,3,1];
[x4,fv4,ef,out4]=fmincon(@pro71,x0,[],[],[],[],v1,v2,@pro71con,opt2)
x0=[1.1,1.1,1.1,1.1];
[x5,fv5,ef,out5]=fmincon(@pro71,x0,[],[],[],[],v1,v2,@pro71con,opt1)
x0=[1.1,1.1,1.1,1.1];
[x6,fv6,ef,out6]=fmincon(@pro71,x0,[],[],[],[],v1,v2,@pro71con,opt2)
[x1;x2;x3;x4;x5;x6]
[fv1;fv2;fv3;fv4;fv5;fv6]
[out1.iterations;out2.iterations;out3.iterations;out4.iterations;out5.iterations;out6.iterations]
[out1.funcCount;out2.funcCount;out3.funcCount;out4.funcCount;out5.funcCount;out6.funcCount]
5. 第一题——（3）小题的非线性约束函数
function[c,ceq,g,geq]=pro71con2(x)
c=[x(1)*x(2)-x(1)-x(2)+1.5;-x(1)*x(2)-10];

ceq=x(1)*x(2)*x(3)*x(4)-16;
ifnargout>2
g=[x(2)-1,-x(2);x(1)-1,-x(1);0,0;0,0];
geq=[x(2)*x(3)*x(4);x(1)*x(3)*x(4);x(1)*x(2)*x(4);x(1)*x(2)*x(3)];
end
end
6. 第一题——（3）小题的求解
v1=[-10,-10,-10,-10];
v2=[10,10,10,10];
A2=[1,1,0,0];
b2=0;
opt1=optimset('largescale','off','MaxIter',3000,'MaxFun',20230);
opt2=optimset(opt1,'GradObj','on','GradConstr','on','DerivativeCheck','on');
x0=[-3,-1,-3,-1];
[x1,fv1,ef,out1]=fmincon(@pro71,x0,[],[],A2,b2,v1,v2,@pro71con2,opt1)
x0=[-3,-1,-3,-1];
[x2,fv2,ef,out2]=fmincon(@pro71,x0,[],[],A2,b2,v1,v2,@pro71con2,opt2)
x0=[3,1,3,1];
[x3,fv3,ef,out3]=fmincon(@pro71,x0,[],[],A2,b2,v1,v2,@pro71con2,opt1)
x0=[3,1,3,1];
[x4,fv4,ef,out4]=fmincon(@pro71,x0,[],[],A2,b2,v1,v2,@pro71con2,opt2)
x0=[1.1,1.1,1.1,1.1];
[x5,fv5,ef,out5]=fmincon(@pro71,x0,[],[],A2,b2,v1,v2,@pro71con2,opt1)
x0=[1.1,1.1,1.1,1.1];
[x6,fv6,ef,out6]=fmincon(@pro71,x0,[],[],A2,b2,v1,v2,@pro71con2,opt2)
[x1;x2;x3;x4;x5;x6]
[fv1;fv2;fv3;fv4;fv5;fv6]
[out1.iterations;out2.iterations;out3.iterations;out4.iterations;out5.iterations;out6.iterations]
[out1.funcCount;out2.funcCount;out3.funcCount;out4.funcCount;out5.funcCount;out6.funcCount]
7. 第二题——临时料场的线性规划
a=[1.258.750.55.7537.25];
b=[1.250.754.7556.57.75];
x=[52];
y=[1,7];
forj=1:2
fori=1:6
s((j-1)*6+i)=sqrt((x(j)-a(i))^2+(y(j)-b(i))^2);
end
end
d=[3547611];
e=[2023];
A1=[;];
b1=e;
A2=zeros(6,12);
fori=1:6
A2(i,i)=1;
A2(i,i+6)=1;
end
b2=d;
v1=zeros(1,12);
opt=optimset('Largescale','off','Simplex','on');
[x,f,exitflag,output,lag]=linprog(s,A1,b1,A2,b2,v1)
x'
8. 第二题——非线性规划的目的函数
functionf=pro72(c,a,b)
forj=1:2
fori=1:6
s((j-1)*6+i)=sqrt((c(11+2*j)-a(i))^2+(c(12+2*j)-b(i))^2);
end
end
f=0;
fori=1:12
f=f+c(i)*s(i);
end
end
9. 第二题——求解非线性规划
a=[1.258.750.55.7537.25];
b=[1.250.7547556.57.75];
d=[3547611];
e=[2023];
A1=[0000;0000];
b1=e;
A2=zeros(6,16);
fori=1:6
A2(i,i)=1;
A2(i,i+6)=1;
end
b2=d;
x0=[zeros(1,12)5127];
v1=zeros(1,16);
opt=optimset('Largescale','off','Maxfun',20230,'Maxiter',3000);
[x,fv,ef,out,lag,grad,hess]=fmincon(@pro72,x0,A1,b1,A2,b2,v1,[],[]opt,a,b)
10. 第二题——消元的目的函数
functionf=pro73(c,a,b,d)
f=0;
fori=1:6
s1=sqrt((c(7)-a(i))^2+(c(8)-b(i))^2);
s2=sqrt((c(9)-a(i))^2+(c(10)-b(i))^2);
f=f+s1*c(i)+s2*(d(i)-c(i));
end
end
11. 第二题——消元的非线性约束的求解
a=[1.258.750.55.7537.25];
b=[1.250.754.7556.57.75];
d=[3547611];
e=20;
A1=[;-1-1-1-1-1-10000];
b1=[ee-sum(d)];
x0=[zeros(1,6)5127];
v1=zeros(1,10);
v2=[d20232023];
opt=optimset('Largescale','off','Maxfun',20230,'Maxiter',3000);
[x,fv,ef,out,lag,grad,hess]=fmincon(@pro73,x0,A1,b1,[],[],v1,v2,[],opt,a,b,d)