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 函数求约束极值问题