数学实验——约束规划

实验7约束优化
分1黄浩2011011743
一、实验目的
1.学会用MATLAB工具箱求解非线性规划的方法。

2.练习建立实际问题的非线性规划模型。

二、实验内容
1.《数学实验》第一版（问题5）
问题叙述：
对问题：
min{100(x2−x12)2+(1−x1)2+90(x4−x32)2+(1−x3)2+ 10.1[(1−x2)2+(1−x4)2]+19.8(x2−1)(x4−1)}，增加以下条件，并分别取初值(-3，-1，-3，-1)和(3，1，3，1)，求解非线性优化：
(1)−10≤x i≤10；
(2)−10≤x i≤10，x1x2−x1−x2+1.5≤0，x1x2+10≥0，−100≤x1x2x3x4≤100；
(3)−10≤x i≤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），所得结果如下：
由上表可见，目标函数的最优值近似为0，最优解近似为(1,1,1,1)，而且不论选取哪个初值、哪种计算方法，最终都能得到近似正确的结果。

另外，因为在命令中没有对误差上限进行设置（即默认值），因此这六个结果的误差数量级都是相同的。

若进行更深入地分析，不难发现，在相同初值和相同精度1条件下，数值方法和分析方法的求解误差和迭代次数相差不大，但分析方法的目标函数调用次数明显少于数值方法。

另一方面，在相同计算方法的条件下，不同初值的迭代次数和调用次数也不同，而且这与初值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），所得结果如下：
因为增加了新的约束条件，最优解和最优值都与上一问有很大出入，而且选取不同的初值也会造成最优解的不同，这是由于matlab只能根据一定的初值和步长求解局部最优点而造成的，获得全局最优点需要枚举所有的局部最优点，并1在上一段中有解释，即精度TolFun和TolX都为相同的默认值
进行比较才能得到。

与上一问相似的地方是，数值方法和分析方法在解的精度和迭代次数上相近，但分析方法的目标函数调用次数较少。

（3）小题
这一小问相对于前者增加了非线性的等式约束和线性的等式约束，因此要修改约束函数的M文件（程序见四.5），分别设初值为(-3,-1,-3,-1)、(3,1,3,1)和(1.1,1.1,1.1,1.1)，使用数值方法和分析方法分别求解（程序见四.6），所得结果如下：
由上表可知，更换约束条件后，题目给定的两个初值(-3,-1,-3,-1)和(3,1,3,1)产生了不同的最优解和最优值，但都能以较少的迭代次数收敛到各自的局部最优解，而且数值方法和分析方法的精度十分接近，只是后者的目标函数调用次数较少。

然而，若使用初值(1.1,1.1,1.1,1.1)，则在目标函数调用次数上限MaxFun=20000的条件下无法完成求解，而且搜索过程是发散的，无法达到最优解。

得出结论：
通过对这个非线性规划问题的求解，我获得了以下几点启示：
1)Matlab数值方法的求解精度很高，一般情况下没有必要使用分析方法求
解梯度，后者只是在目标函数的调用次数上较少，对于精度的贡献比较
小，在中小规模运算和迭代次数较少时可直接使用数值方法。

2)在最优解唯一的情况下，使用不同初值（要求收敛）所得最优解和最优
值近似，但迭代次数有较大差异，而且迭代次数与初值x0-最优解x之
间的距离没有本质联系。

而在最优解不唯一的情况下，不同初值可能会
产生不同的最优解，它们都是局部极值点，需要比较所有局部极值点后
才能得到全局极值点（最值点）。

3)在优化实验中，对优化函数的控制参数进行设置看似多余，实际上很有
必要，例如设置最大调用次数和最大迭代次数来防止发散情况下的无限
迭代，以及调节输出精度来满足更高的要求等等。

4)这次进行的是小规模的非线性优化实验，而在更深层次的工程领域，可
能会出现大规模的实验，求解速度会十分缓慢，这时，使用更好的算法
和更合适的初值就变得尤为重要。

2.《数学实验》第一版（例题：供应与选址）
问题描述：
某公司有6个建筑工地要开工，每个工地的位置(用平面坐标x，y表示，距离单位：km)及水泥日用量d(t)由下表给出。

目前有两个临时料场位于
A(5，1)，B(2，7)，日储量各有20t，假设从料场到工地之间均有直线道路相连，试制订每天的供应计划，即从两料场分别向各工地运送多少吨水泥，使总的吨公里数最小。

为进一步减少吨公里数，打算舍弃两个临时料场，改建两个新的，日储量仍各为20吨，问应建在何处？节省的吨公里数有多大？
工地的位置(x，y)及水泥日用量d
模型转换及实验过程：
设：各工地的位置分别为(a i,b i)，水泥的日用量分别为d i（i=1,…,6，分别表示六个工地），料场位置设为(x j,y j)，日储量分别为e j（j=1,2，分别表示A,B料场），从料场j运往工地i的水泥运送量为c ij，则原问题可以表示为以下的优化问题：
min f=∑∑c ij
6
i=1√(x
j
−a i)2+(y j−b i)2
2
j=1
s.t.∑c ij≤e j,j=1,2
6
i=1
2
∑c ij
=d i,i=1,2,…,6
j=1
c ij≥0,i=1,2,……,6,j=1,2
以c ij作为决策变量，不难发现这是一个线性规划问题，使用matlab求解这个线性规划（程序见四.7），结果如下：
四舍五入后得到：
因此，当各料场运送量按上表执行时，总的吨公里数为136.228，且为全局最小值。

如果舍弃现有的A、B料场，改建新料场，则现在的优化问题变成了非线性优化问题，目标函数变成了非线性的，决策变量是料场向工地的运送吨数和料场位置，而约束条件及上下界都不变。

先编写目标函数的M文件（程序见四.8），将原A、B料场的位置作为新优化问题的初值，然后使用matlab求解这个非线性约束（程序见四.9），结果如下（经过了四舍五入）：
在这种情况的下的总吨公里数为88.883，迭代次数为120次。

由于matlab的非线性优化只能通过一个初值来求得局部极小值，因此我们换用其他的位置初值进行求解（运送量初值均设为0,程序略），结果如下：
由上表可见，选取不同的初值都能够收敛到一定的局部最优解，然而它们的目标函数最优值却各不相同，在当前所有的局部最优解中，最小的为85.266，即上表中的第一个初值所得结果。

与此同时，我们发现∑c ij 2j=1=d i ,i =1,2,…,6是六个线性的等式约束，因此一定可以通过消元法消除六个变量，恰好这六个等式约束都为二元变量，消元就更加容易了。

重新对这个非线性规划问题进行整理，得：
min ∑{c i √(x 1−a i )2+(y 1−b i )2+(d i −c i )√(x 2−a i )2+(y 2−b i )26
i=1
}
s.t.0≤c i ≤d i ,i =1,2,……,6
∑c i ≤e 16
i=1
∑(d i −c i )≤e 26
i=1
因此，整理后的决策变量为新料场A’向六个工地的运送量c i 和两个新料场的位置(x j ,y j )，重新编译目标函数的M 文件（程序见四.10），使用原料场的位置为位置初值，然后用matlab 求解这个非线性规划（程序见四.11），结果如下：
由上表可见，在重新整理了这个优化模型之后，得到了更优化的解，而且发现这种情况下的迭代次数为36次，远小于一开始的120次，此外，通过与前面获得的局部最优值比较可得，85.266是当前获得的全局最优值。

这说明，降低问题的维数对于提高解精度以及减小迭代次数方面有好处的。

当然，这也要视情况而定，当等式约束皆为非线性约束、或是多元的线性约束时，盲目进行消元降维会增加人的计算复杂度，编程时也容易出错，因此，我认为，只有像本题这种线性的二元约束，才有消元的必要，其他情况下盲目消元是无益的，毕竟matlab 的数值功能还是很强大的。

得出结论：
通过实验可得，若使用临时料场A 、B ，则该问题为线性规划问题，当料场
向各个工地的运送量如下表所示时，总的吨公里数最少：
若改建新料场，则当料场位置和运送量如下表所示时，总的吨公里数最少：
而且通过这两个表的比较可见，新料场的吨公里数比临时料场的吨公里数减少了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);
if nargout>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',20000);
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.fun cCount;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];
if nargout>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) *x(4),-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',20000);
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];
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.fun cCount;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;
if nargout>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',20000);
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];
[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.fun cCount;out6.funcCount]
7.第二题——临时料场的线性规划
a=[1.258.750.55.7537.25];
b=[1.250.754.7556.57.75];
x=[52];
y=[1,7];
for j=1:2
for i=1:6
s((j-1)*6+i)=sqrt((x(j)-a(i))^2+(y(j)-b(i))^2);
end
end
d=[3547611];
e=[2020];
A1=[111111000000;000000111111];
b1=e;
A2=zeros(6,12);
for i=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.第二题——非线性规划的目标函数
function f=pro72(c,a,b)
for j=1:2
for i=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;
for i=1:12
f=f+c(i)*s(i);
end
end
9.第二题——求解非线性规划
a=[1.258.750.55.7537.25];
b=[1.250.754.7556.57.75];
d=[3547611];
e=[2020];
A1=[1111110000000000;0000001111110000];
b1=e;
A2=zeros(6,16);
for i=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',20000,'Maxiter',3000);
[x,fv,ef,out,lag,grad,hess]=fmincon(@pro72,x0,A1,b1,A2,b2,v1,[],[],op t,a,b)
10.第二题——消元的目标函数
function f=pro73(c,a,b,d)
f=0;
for i=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=[1111110000;-1-1-1-1-1-10000];
b1=[ee-sum(d)];
x0=[zeros(1,6)5127];
v1=zeros(1,10);
v2=[d2*******];
opt=optimset('Largescale','off','Maxfun',20000,'Maxiter',3000);
[x,fv,ef,out,lag,grad,hess]=fmincon(@pro73,x0,A1,b1,[],[],v1,v2,[],op t,a,b,d)。