Algorithm之PrA：PrA之nLP整数规划算法经典案例剖析+Matlab编程实现

目录

有约束非线性规划案例分析

1、投资决策问题

2、利用Matlab实现求解下列非线性规划​

无约束极值问题案例分析

1、解析法中的梯度法

2、解析法中的牛顿法

3、Matlab 求无约束极值问题

4、求多元函数的极值

5、罚函数法求解非线性规划

二次规划案例分析

1、求解二次规划

Matlab求约束极值问题

1、fminbnd 函数求约束极值问题​

2、fseminf 函数求约束极值问题​

3、fminimax 函数求约束极值问题​

4、利用梯度求解约束优化问题

有约束非线性规划案例分析

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];
[-meta">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;
[-meta">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-emphasis">*x(2)^4+x(1)^2*x(2)^2;
df=[4-emphasis">*x(1)^3+2*x(1)-emphasis">*x(2)^2;100*x(2)^3+2-emphasis">*x(1)^2*x(2)];
d2f=[2-emphasis">*x(1)^2+2*x(2)^2,4-emphasis">*x(1)*x(2)
-code">    4*x(1)*x(2),300*x(2)^2+2*x(1)^2];

(3)、编写主程序文件example5.m 如下

clc
x=[2;2];
[-meta">f0,g1,g2]=nwfun(x);
while norm(g1)>0.00001
    p=-inv(g2)*g1;
    x=x+p;
    [-meta">f0,g1,g2]=nwfun(x);
end
x, f0

(4)、扩展：如果目标函数是非二次函数，一般地说，用Newton 法通过有限轮迭代并不能保证可求得其最优解。为了提高计算精度，我们在迭代时可以采用变步长计算上述问题，编写主程序文件

example5_2 如下：
clc,clear
x=[2;2];
[-meta">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;
[-meta">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(-variable">@fun4,x0)

5、罚函数法求解非线性规划

(1)、编写 M 文件 test.m

function g-operator">=test(x);
M-operator">=50000;
f-operator">=x(1)-operator">^2-operator">+x(2)-operator">^2-operator">+8;
g-operator">=f-operator">-M-operator">*min(x(1),0)-operator">-M-operator">*min(x(2),0)-operator">-M-operator">*min(x(1)-operator">^2-operator">-x(2),0)-operator">+-operator">...
    M-operator">*abs(-operator">-x(1)-operator">-x(2)-operator">^2-operator">+2);
 
 
 
或者是利用Matlab的求矩阵的极小值和极大值函数编写test.m如下：
function g-operator">=test(x);
M-operator">=50000;
f-operator">=x(1)-operator">^2-operator">+x(2)-operator">^2-operator">+8;
g-operator">=f-operator">-M-operator">*sum(min([x';zeros(1,2)]))-operator">-M-operator">*min(x(1)-operator">^2-operator">-x(2),0)-operator">+-operator">...
    M-operator">*abs(-operator">-x(1)-operator">-x(2)-operator">^2-operator">+2);
 
 
我们也可以修改罚函数的定义，编写test.m如下：
function g-operator">=test(x);
M-operator">=50000;
f-operator">=x(1)-operator">^2-operator">+x(2)-operator">^2-operator">+8;
g-operator">=f-operator">-M-operator">*min(min(x),0)-operator">-M-operator">*min(x(1)-operator">^2-operator">-x(2),0)-operator">+M-operator">*(-operator">-x(1)-operator">-x(2)-operator">^2-operator">+2)-operator">^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];
[-meta">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))-emphasis">*(4*x(1)^2+2-emphasis">*x(2)^2+4*x(1)-emphasis">*x(2)+2*x(2)+1);
df=[exp(x(1))-emphasis">*(4*x(1)^2+2-emphasis">*x(2)^2+4*x(1)-emphasis">*x(2)+8*x(1)+6-emphasis">*x(2)+1);exp(x(1))*(4-emphasis">-emphasis">*x(2)
-emphasis">+4*x(1)+2)];

(2)、编写M 文件fun10.m 定义约束条件及约束条件的梯度函数：

function -punctuation">[c-punctuation">,ceq-punctuation">,dc-punctuation">,dceq-punctuation">]-operator">=fun10-punctuation">(x-punctuation">);
c-operator">=-punctuation">[x-punctuation">(1-punctuation">)-operator">*x-punctuation">(2-punctuation">)-operator">-x-punctuation">(1-punctuation">)-operator">-x-punctuation">(2-punctuation">)-operator">+1.5;-operator">-x-punctuation">(1-punctuation">)-operator">*x-punctuation">(2-punctuation">)-operator">-10-punctuation">];
dc-operator">=-punctuation">[x-punctuation">(2-punctuation">)-operator">-1-punctuation">,-operator">-x-punctuation">(2-punctuation">);x-punctuation">(1-punctuation">)-operator">-1-punctuation">,-operator">-x-punctuation">(1-punctuation">)-punctuation">];
ceq-operator">=-punctuation">[-punctuation">];dceq-operator">=-punctuation">[-punctuation">];

(3)、调用函数fmincon，编写主函数文件example13.m 如下

%采用标准算法
options=optimset('largescale','off');
%采用梯度
options=optimset(options,'GradObj','on','GradConstr','on');
[x,y]=fmincon(@fun9,rand(2,1),[],[],[],[],[],[],@fun10,options)