Algorithm之PrA:PrA之nLP非线性规划算法经典案例剖析+Matlab编程实现
Algorithm之PrA:PrA之nLP整数规划算法经典案例剖析+Matlab编程实现
有约束非线性规划案例分析
1、投资决策问题
某企业有n 个项目可供选择投资,并且至少要对其中一个项目投资。已知该企业拥有总资金 A元,投资于第i(i = 1…,n)个项目需花资金 ai元,并预计可收益i b 元。试选择最佳投资方案。
(1)、根据已知列出数学公式
(2)、分析问题,列出问题模型
(3)、分析模型,寻找求解模型方法
2、利用Matlab实现求解下列非线性规划
(1)、编写M 文件fun1.m 定义目标函数
function f=fun1(x);
f=sum(x.^2)+8;
(2)、编写M文件fun2.m定义非线性约束条件
function [g,h]=fun2(x);
g=[-x(1)^2+x(2)-x(3)^2
x(1)+x(2)^2+x(3)^3-20]; %非线性不等式约束
h=[-x(1)-x(2)^2+2
x(2)+2*x(3)^2-3]; %非线性等式约束
(3)、编写主程序文件example2.m 如下
options=optimset('largescale','off');
[x,y]=fmincon('fun1',rand(3,1),[],[],[],[],zeros(3,1),[], ...
'fun2', options)
(4)、可以求得
无约束极值问题案例分析
1、解析法中的梯度法
用最速下降法求解无约束非线性规划问题
(1)、求导
编写 M 文件detaf.m,定义函数 f (x)及其梯度列向量如下
function [f,df]=detaf(x);
f=x(1)^2+25*x(2)^2;
df=[2*x(1)
50*x(2)];
(2)、编写主程序文件zuisu.m如下
clc
x=[2;2];
[f0,g]=detaf(x);
while norm(g)>0.000001
p=-g/norm(g);
t=1.0;f=detaf(x+t*p);
while f>f0
t=t/2;
f=detaf(x+t*p);
end
x=x+t*p;
[f0,g]=detaf(x);
end
x,f0
2、解析法中的牛顿法
用Newton 法求解
(1)、求梯度
(2)、编写 M 文件nwfun.m 如下:
编写 M 文件nwfun.m 如下:
function [f,df,d2f]=nwfun(x);
f=x(1)^4+25*x(2)^4+x(1)^2*x(2)^2;
df=[4*x(1)^3+2*x(1)*x(2)^2;100*x(2)^3+2*x(1)^2*x(2)];
d2f=[2*x(1)^2+2*x(2)^2,4*x(1)*x(2)
4*x(1)*x(2),300*x(2)^2+2*x(1)^2];
(3)、编写主程序文件example5.m 如下
clc
x=[2;2];
[f0,g1,g2]=nwfun(x);
while norm(g1)>0.00001
p=-inv(g2)*g1;
x=x+p;
[f0,g1,g2]=nwfun(x);
end
x, f0
(4)、扩展:如果目标函数是非二次函数,一般地说,用Newton 法通过有限轮迭代并不能保证可求得其最优解。为了提高计算精度,我们在迭代时可以采用变步长计算上述问题,编写主程序文件
example5_2 如下:
clc,clear
x=[2;2];
[f0,g1,g2]=nwfun(x);
while norm(g1)>0.00001
p=-inv(g2)*g1;p=p/norm(p);
t=1.0;f=nwfun(x+t*p);
while f>f0
t=t/2;f=nwfun(x+t*p);
end
x=x+t*p;
[f0,g1,g2]=nwfun(x);
end
x,f0
3、Matlab 求无约束极值问题
(1)、编写M 文件fun2.m 如下
function [f,g]=fun2(x);
f=100*(x(2)-x(1)^2)^2+(1-x(1))^2;
g=[-400*x(1)*(x(2)-x(1)^2)-2*(1-x(1));200*(x(2)-x(1)^2)];
(2)、编写主函数文件example6.m如下
options = optimset('GradObj','on');
[x,y]=fminunc('fun2',rand(1,2),options)
即可求得函数的极小值。
在求极值时,也可以利用二阶导数,编写 M 文件fun3.m 如下:
function [f,df,d2f]=fun3(x);
f=100*(x(2)-x(1)^2)^2+(1-x(1))^2;
df=[-400*x(1)*(x(2)-x(1)^2)-2*(1-x(1));200*(x(2)-x(1)^2)];
d2f=[-400*x(2)+1200*x(1)^2+2,-400*x(1)
-400*x(1),200];
(3)、编写主函数文件example62.m如下:即可求得函数的极小值。
options = optimset('GradObj','on','Hessian','on');
[x,y]=fminunc('fun3',rand(1,2),options)
4、求多元函数的极值
(1)、编写 f (x)的 M 文件 fun4.m如下
function f=fun4(x);
f=sin(x)+3;
(2)、编写主函数文件example7.m如下,即求得在初值2 附近的极小点及极小值。
x0=2;
[x,y]=fminsearch(@fun4,x0)
5、罚函数法求解非线性规划
(1)、编写 M 文件 test.m
function g=test(x);
M=50000;
f=x(1)^2+x(2)^2+8;
g=f-M*min(x(1),0)-M*min(x(2),0)-M*min(x(1)^2-x(2),0)+...
M*abs(-x(1)-x(2)^2+2);
或者是利用Matlab的求矩阵的极小值和极大值函数编写test.m如下:
function g=test(x);
M=50000;
f=x(1)^2+x(2)^2+8;
g=f-M*sum(min([x';zeros(1,2)]))-M*min(x(1)^2-x(2),0)+...
M*abs(-x(1)-x(2)^2+2);
我们也可以修改罚函数的定义,编写test.m如下:
function g=test(x);
M=50000;
f=x(1)^2+x(2)^2+8;
g=f-M*min(min(x),0)-M*min(x(1)^2-x(2),0)+M*(-x(1)-x(2)^2+2)^2;
(2)、在Matlab 命令窗口输入 [x,y]=fminunc('test',rand(2,1)) 即可求得问题的解。
二次规划案例分析
1、求解二次规划
(1)、编写程序
h=[4,-4;-4,8];
f=[-6;-3];
a=[1,1;4,1];
b=[3;9];
[x,value]=quadprog(h,f,a,b,[],[],zeros(2,1))
Matlab求约束极值问题
1、fminbnd 函数求约束极值问题
解:编写 M 文件fun5.m
function f=fun5(x);
f=(x-3)^2-1;
在Matlab 的命令窗口输入[x,y]=fminbnd('fun5',0,5) 即可求得极小点和极小值。
2、fseminf 函数求约束极值问题
(1)、编写M 文件fun6.m 定义目标函数如下
function f=fun6(x,s);
f=sum((x-0.5).^2);
(2)、编写M 文件fun7.m 定义约束条件如下
function [c,ceq,k1,k2,s]=fun7(x,s);
c=[];ceq=[];
if isnan(s(1,1))
s=[0.2,0;0.2 0];
end
%取样值
w1=1:s(1,1):100;
w2=1:s(2,1):100;
%半无穷约束
k1=sin(w1*x(1)).*cos(w1*x(2))-1/1000*(w1-50).^2-sin(w1*x(3))-x(3)-1;
k2=sin(w2*x(2)).*cos(w2*x(1))-1/1000*(w2-50).^2-sin(w2*x(3))-x(3)-1;
%画出半无穷约束的图形
plot(w1,k1,'-',w2,k2,'+');
(3)调用函数fseminf,在Matlab 的命令窗口输入即可。
[x,y]=fseminf(@fun6,rand(3,1),2,@fun7)
3、fminimax 函数求约束极值问题
(1)、编写M 文件fun8.m 定义向量函数如下
function f=fun8(x);
f=[2*x(1)^2+x(2)^2-48*x(1)-40*x(2)+304
-x(1)^2-3*x(2)^2
x(1)+3*x(2)-18
-x(1)-x(2)
x(1)+x(2)-8];
(2)、调用函数fminimax
[x,y]=fminimax(@fun8,rand(2,1))
4、利用梯度求解约束优化问题
分析:当使用梯度求解上述问题时,效率更高并且结果更准确。题目中目标函数的梯度为
(1)、编写M 文件fun9.m 定义目标函数及梯度函数
function [f,df]=fun9(x);
f=exp(x(1))*(4*x(1)^2+2*x(2)^2+4*x(1)*x(2)+2*x(2)+1);
df=[exp(x(1))*(4*x(1)^2+2*x(2)^2+4*x(1)*x(2)+8*x(1)+6*x(2)+1);exp(x(1))*(4*x(2)
+4*x(1)+2)];
(2)、编写M 文件fun10.m 定义约束条件及约束条件的梯度函数:
function [c,ceq,dc,dceq]=fun10(x);
c=[x(1)*x(2)-x(1)-x(2)+1.5;-x(1)*x(2)-10];
dc=[x(2)-1,-x(2);x(1)-1,-x(1)];
ceq=[];dceq=[];
(3)、调用函数fmincon,编写主函数文件example13.m 如下
%采用标准算法
options=optimset('largescale','off');
%采用梯度
options=optimset(options,'GradObj','on','GradConstr','on');
[x,y]=fmincon(@fun9,rand(2,1),[],[],[],[],[],[],@fun10,options)