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];
[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)

 

 

 

线性规划问题
一、线性规划问题及其标准型
1.问题的求解过程







2.标准形式




二、非标准形问题的转换

三、线性规划问题中的Matlab求解



四、案例分析


问题分析



分析结果







整理出表达式



代码实现

可以把下面代码，先存到M文件，然后运行，也可以逐行粘贴至命令行运行
f=[-2.85,-3.05,-2.9,3.1,3.25,2.95];
f=f';f=-f;%价值向量，并转化为求最小值时对应的价值向量
A=[0,0,0,1,0,0;
  -1,0,0,1,1,0;
  -1,-1,0,1,1,1;
  1,0,0,-1,0,0;
 1,1,0,-1,-1,0;
 2.85,0,0,-3.1,0,0;
 2.85,3.05,0,-3.1,-3.25,0;
 2.85,3.05,2.9,-3.1,-3.25,-2.95];
b=[1000;1000;1000;4000;4000;20000;20000;20000];%线性不等式约束
aeq=[1,1,1,-1,-1,-1];beq=1000;%线性等式约束
lb=zeros(6,1);%函数值下限
[x,y]=linprog(f,A,b,aeq,beq,lb);%输出x中后三项分别代表y1,y2,y3
 x=reshape(x,1,6),-y%重新组合矩阵x，取负求得原本的最大值

运行结果如下：


五、蒙特卡洛相关方法的介绍

x=unifrnd(0,12,[1,10000000]);
y=unifrnd(0,9,[1,10000000]);
Pinshu=sum(y<x.^2&x<=3)+sum(y<12-x&x>=3);
Area_appr=12*9*Pinshu/10^7

%运行结果在49.5附近，由于是随机模拟，因此每次的结果都是不一样的。

Algorithm之PrA：PrA之LP线性规划算法经典案例剖析

 

 

目录

一、以例题分析步骤理解LP线性规划算法

二、Matlab编程实现

三、将问题可以转化为线性规划的问题

四、建模思路全过程

 

 

 

一、以例题分析步骤理解LP线性规划算法

1、分析问题

      某机床厂生产甲、乙两种机床，每台销售后的利润分别为4000 元与3000 元。生产甲机床需用 A、B机器加工，加工时间分别为每台 2 小时和 1 小时；生产乙机床需用 A、B、C三种机器加工，加工时间为每台各一小时。若每天可用于加工的机器时数分别为A 机器10 小时、B 机器8 小时和C 机器7 小时，问该厂应生产甲、乙机床各几台，才能使总利润最大？

分析：上述问题的数学模型：设该厂生产x1 台甲机床和x2 台乙机床时总利润最大，则x1、x2应满足


      这里变量x1、x2称之为决策变量，（1）式被称为问题的目标函数，（2）中的几个不等式是问题的约束条件，记为s.t.(即subject to)。由于上面的目标函数及约束条件均为线性函数，故被称为线性规划问题。
经验：总之，线性规划问题是在一组线性约束条件的限制下，求一线性目标函数最大或最小的问题。在解决实际问题时，把问题归结成一个线性规划数学模型是很重要的一步，但往往也是困难的一步，模型建立得是否恰当，直接影响到求解。而选适当的决策变量，是我们建立有效模型的关键之一。

2、线性规划的 Matlab 标准形式

         线性规划的目标函数可以是求最大值，也可以是求最小值，约束条件的不等号可以是小于号也可以是大于号。为了避免这种形式多样性带来的不便，Matlab 中规定线性规划的标准形式为

                  

    其中 c 和x 为n 维列向量， A 、Aeq 为适当维数的矩阵， b 、beq 为适当维数的列向量。

3 线性规划问题的解的概念

      一般线性规划问题的（数学）标准型为



可行解 满足约束条件（4）的解x=(x1,x2,x3……xn)，称为线性规划问题的可行解，而使目标函数（3）达到最大值的可行解叫最优解。
可行域 所有可行解构成的集合称为问题的可行域，记为R 。

4 线性规划的图解法

        图解法简单直观，有助于了解线性规划问题求解的基本原理。我们先应用图解法来求解例1。对于每一固定的值z ，使目标函数值等于z 的点构成的直线称为目标函数等位线，当z 变动时，我们得到一族平行直线。对于例1，显然等位线越趋于右上方，其上的点具有越大的目标函数值。不难看出，本例的最优解为x* = (2,6)T ，最优目标值z* = 26。

5 求解线性规划的Matlab 解法

        单纯形法是求解线性规划问题的最常用、最有效的算法之一。这里我们就不介绍单纯形法，有兴趣的读者可以参看其它线性规划书籍。下面我们介绍线性规划的Matlab解法。


      基本函数形式为 linprog(c,A,b)，它的返回值是向量x 的值。还有其它的一些函数调用形式（在 Matlab 指令窗运行 help linprog 可以看到所有的函数调用形式），如：


     这里fval 返回目标函数的值，LB 和UB 分别是变量x 的下界和上界， x0是x 的初始值，OPTIONS 是控制参数。

二、Matlab编程实现

1、求解下列线性规划问题01


解：编写M 文件 → 将M文件存盘，并命名为example1.m → 在Matlab指令窗运行example1即可得所求结果。

c=[2;3;-5];
a=[-2,5,-1;1,3,1]; b=[-10;12];
aeq=[1,1,1];
beq=7;
x=linprog(-c,a,b,aeq,beq,zeros(3,1))
value=c'*x


2、求解线性规划问题02



解：编写Matlab程序如下：

c=[2;3;1];
a=[1,4,2;3,2,0];
b=[8;6];
[x,y]=linprog(c,-a,-b,[],[],zeros(3,1))

三、将问题可以转化为线性规划的问题

很多看起来不是线性规划的问题也可以通过变换变成线性规划的问题来解决

3.1、比如将带有绝对值的复杂问题，转化为LP问题




四、建模思路全过程

1、提出问题：投资的收益和风险问题





2、符号规定和基本假设



3、模型的分析与建立



4、模型一的求解



clc,clear
a=0;
hold on
while a<0.05
c=[-0.05,-0.27,-0.19,-0.185,-0.185];
A=[zeros(4,1),diag([0.025,0.015,0.055,0.026])];
b=a*ones(4,1);
Aeq=[1,1.01,1.02,1.045,1.065];
beq=1;
LB=zeros(5,1);
[x,Q]=linprog(c,A,b,Aeq,beq,LB);
Q=-Q;
plot(a,Q,'*r');
a=a+0.001;
end
xlabel('a'),ylabel('Q')

5、结果分析



 

 

 

 

 

相关文章参考

《数学建模算法与程序大全》

 

 

 

文章目录
题目
问题简化
约束条件
目标函数
主函数
目标函数
约束条件
fmincon函数详解

题目
华安机械厂的厂长正考虑将该厂的一部分在市区的生产车间搬至该市的卫星城镇，好处是土地、房租费及排污处理费用等都较便宜，但这样做会增加车间之间的交通运输费用。

该厂原在市区车间有A、B、C、D、E 五个，计划搬迁去的卫星城镇有甲、乙两处。规定无论留在市区或甲、乙两卫星城镇均不得多于3 个车间。



但搬迁后带来运输费用增加由$C_{ik}$ 和$D_{jk}$ 值决定， $C_{ik}$为$i$和$k$车间之间的年运量， $D_{jk}$ 为市区同卫星城镇

间单位运量的运费，具体数据分别见下表



请为厂长提供一个决策建议方案，哪几个车间搬至卫星城镇及搬至甲还是乙，能带来最大的经济上的好处。
问题简化
该问题其实本质上为：非线性0-1规划问题，可以参考小黄书中的0-1整数规划和非线性规划内容。

首先，设置决策变量$X_{ij}$为第$i$个车间的第$j$种状态，其中 $i = 1，2，3，4，5$分别表示车间$A、B、C、D、E$ ，而$j$的$i = 1，2，3$表示为搬去甲地，搬乙地，留在市区，且xij 为0,1 变量，表示不搬、搬。如下所示：






搬甲地
搬乙地
留在市区




A
X11
X12
X13


B
X21
X22
X23


C
X31
X32
X33


D
X41
X42
X43


E
X51
X52
X53


约束条件
各车间最终处于3 种状态之一，有$\sum_{j=1}^3X_{ij}=1,i=1,2,3,4,5$

即：

$X11+X12+X13=1;$

$X21+X22+X23=1;$

$X31+X32+X33=1;$

$X41+X42+X43=1;$

$X51+X52+X53=1;$

同理，又由于规定无论留在市区或甲、乙两卫星城镇均不得多于3 个车间，有$\sum_{i=1}^5X_{ij}=1,j=1,2,3$

即：
$X11+X21+X31+X41+X51<=3;$

$X12+X22+X32+X42+X52<=3;$

$X13+X23+X33+X43+X53<=3;$
目标函数
最大经济效益 ：$节约费用*10000-运输费用$

节约费用：$100*x11+100*x12+150*x21+200*x22+100*x31+150*x32+200*x41+150*x42+50*x51+150*x52$

运输费用：$1000*b13+1500*b14+1400*b23+1200*b24+2000*b35+700*b45$
为了得到最大经济效益，需要使节约费用最大以及搬迁后的运输费用最小。故可建立相应模型，模型这里结合表格以及上述约束条件便可，我这里就不再写一次了。下面是matlab代码：
主函数
X0 = zeros(15, 1);
A = [];
B = [];
Aeq = [];
Beq = [];
LB = zeros(15, 1);
UB = ones(15,1);

[x, fval] = fmincon('objfun', X0, A, B, Aeq, Beq, LB, UB, 'NONLCON')   % fmincon函数文后有详解。

目标函数
function f=objfun(x)
f=-10000*(100*x(1)+100*x(2)+150*x(4)+200*x(5)+100*x(7)+150*x(8)+200*x(10)+150*x(11)+50*x(13)+150*x(14))+...
    (1000*(50*x(1)*x(7)+140*x(1)*x(8)+130*x(1)*x(9)+50*x(2)*x(8)+90*x(2)*x(9)+50*x(3)*x(9))+...
    1500*(50*x(1)*x(10)+140*x(1)*x(11)+130*x(1)*x(12)+50*x(2)*x(11)+90*x(2)*x(12)+50*x(3)*x(12))+...
    1400*(50*x(4)*x(7)+140*x(4)*x(8)+130*x(4)*x(9)+50*x(5)*x(8)+90*x(5)*x(9)+50*x(6)*x(9))+...
    1200*(50*x(4)*x(10)+140*x(4)*x(11)+130*x(4)*x(12)+50*x(5)*x(11)+90*x(5)*x(12)+50*x(6)*x(12))+...
    2000*(50*x(7)*x(13)+140*x(7)*x(14)+130*x(7)*x(15)+50*x(8)*x(14)+90*x(8)*x(15)+50*x(9)*x(15))+...
    700*(50*x(10)*x(13)+140*x(10)*x(14)+130*x(10)*x(15)+50*x(11)*x(14)+90*x(11)*x(15)+50*x(12)*x(15)));
end

约束条件
function [g, ceq] = NONLCON(x)
 g = [x(1)+x(4)+x(7)+x(10)+x(13)-3
     x(2)+x(5)+x(8)+x(11)+x(14)-3
     x(3)+x(6)+x(9)+x(12)+x(15)-3];   %非线性不等式约束
 ceq = [x(1)+x(2)+x(3)-1
        x(4)+x(5)+x(6)-1
        x(7)+x(8)+x(9)-1
        x(10)+x(11)+x(12)-1
        x(13)+x(14)+x(15)-1];        %非线性等式约束
end

fmincon函数详解

在matlab中，fmincon函数可以求解带约束的非线性多变量函数(Constrained nonlinear multivariable function)的最小值,即可以用来求解非线性规划问题

matlab中，非线性规划模型的写法如下

[x,fval]=fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon,options)


x的返回值是决策向量x的取值，fval的返回值是目标函数f(x)的取值


fun是用M文件定义的函数f(x),代表了(非)线性目标函数


x0是x的初始值


A,b,Aeq,beq定义了线性约束 ,如果没有线性约束，则A=[],b=[],Aeq=[],beq=[]


lb和ub是变量x的下界和上界，如果下界和上界没有约束，则lb=[],ub=[],也可以写成lb的各分量都为 -inf,ub的各分量都为inf


nonlcon是用M文件定义的非线性向量函数约束


options定义了优化参数，不填写表示使用Matlab默认的参数设置