MATLAB优化应用非线性规划.doc

MATLAB优化应用
§1 线性规划模型
一、线性规划课题：
实例1：生产计划问题
假设某厂计划生产甲、乙两种产品，现库存主要材料有A类3600公斤，B类2000公斤，C类3000公斤。每件甲产品需用材料A类9公斤，B类4公斤，C类3公斤。每件乙产品，需用材料A类4公斤，B类5公斤，C类10公斤。甲单位产品的利润70元，乙单位产品的利润120元。问如何安排生产，才能使该厂所获的利润最大。
建立数学模型：
设x1、x2分别为生产甲、乙产品的件数。f为该厂所获总润。
       max f=70x1+120x2
       s.t  9x1+4x2≤3600
        4x1+5x2≤2000
        3x1+10x2≤3000
        x1,x2≥0
实例2：投资问题
某公司有一批资金用于4个工程项目的投资，其投资各项目时所得的净收益(投入资金锪百分比)如下表：
工程项目收益表
工程项目
A
B
C
D
收益(%)
15
10
8
12
由于某种原因，决定用于项目A的投资不大于其他各项投资之和而用于项目B和C的投资要大于项目D的投资。试确定全文该公司收益最大的投资分配方案。
建立数学模型：
设x1、 x2 、x3 、x4分别代表用于项目A、B、C、D的投资百分数。
      max f=0.15x1+0.1x2+0.08 x3+0.12 x4
       s.t  x1-x2- x3- x4≤0
        x2+ x3- x4≥0
              x1+x2+x3+ x4=1
              xj≥0  j=1,2,3,4
实例3：运输问题
有A、B、C三个食品加工厂，负责供给甲、乙、丙、丁四个市场。三个厂每天生产食品箱数上限如下表：
工厂
A
B
C
生产数
60
40
50
四个市场每天的需求量如下表：
市场
甲
乙
丙
丁
需求量
20
35
33
34
从各厂运到各市场的运输费(元/每箱)由下表给出：
收

点

发

点

市  场
甲
乙
丙
丁
工
厂
A
2
1
3
2
B
1
3
2
1
C
3
4
1
1
求在基本满足供需平衡的约束条件下使总运输费用最小。
建立数学模型：
设ai j为由工厂i运到市场j的费用，xi j 是由工厂i运到市场j的箱数。bi是工厂i的产量，dj是市场j的需求量。

b= ( 60 40 50 )   d= ( 20 35 33 34 )

       s.t 

        x i j≥0

当我们用MATLAB软件作优化问题时，所有求maxf 的问题化为求min(-f )来作。约束g i (x)≥0，化为 -g i≤0来作。
上述实例去掉实际背景，归结出规划问题：目标函数和约束条件都是变量x的线性函数。
形如：    (1)        min f T X
                     s.t  A X≤b
       Aeq X =beq
lb≤X≤ub
    其中X为n维未知向量，f T=[f1,f2,…fn]为目标函数系数向量，小于等于约束系数矩阵A为m×n矩阵，b为其右端m维列向量，Aeq为等式约束系数矩阵，beq为等式约束右端常数列向量。lb,ub为自变量取值上界与下界约束的n维常数向量。
二．线性规划问题求最优解函数：
       调用格式：  x=linprog(f,A,b)
                            x=linprog(f,A,b,Aeq,beq)
                            x=linprog(f,A,b,Aeq,beq,lb,ub)
                            x=linprog(f,A,b,Aeq,beq,lb,ub,x0)
                            x=linprog(f,A,b,Aeq,beq,lb,ub,x0,options)
              [x,fval]=linprog(…)
              [x, fval, exitflag]=linprog(…)
              [x, fval, exitflag, output]=linprog(…)
              [x, fval, exitflag, output, lambda]=linprog(…)
       说明：x=linprog(f,A,b)返回值x为最优解向量。
       x=linprog(f,A,b,Aeq,beq) 作有等式约束的问题。若没有不等式约束，则令A=[ ]、b=[ ] 。
       x=linprog(f,A,b,Aeq,beq,lb,ub,x0,options) 中lb ,ub为变量x的下界和上界，x0为初值点，options为指定优化参数进行最小化。
Options的参数描述：
Display   显示水平。 选择’off’ 不显示输出；选择’iter’显示每一 步迭代过程的输出；选择’final’ 显示最终结果。
MaxFunEvals 函数评价的最大允许次数
Maxiter 最大允许迭代次数
TolX   x处的终止容限     
       [x,fval]=linprog(…) 左端 fval 返回解x处的目标函数值。
[x,fval,exitflag,output,lambda]=linprog(f,A,b, Aeq,beq,lb,ub,x0) 的输出部分：
exitflag 描述函数计算的退出条件：若为正值，表示目标函数收敛于解x处；若为负值，表示目标函数不收敛；若为零值，表示已经达到函数评价或迭代的最大次数。
output 返回优化信息：output.iterations表示迭代次数；output.algorithm表示所采用的算法；outprt.funcCount表示函数评价次数。
lambda 返回x处的拉格朗日乘子。它有以下属性：
       lambda.lower-lambda的下界；
       lambda.upper-lambda的上界；
       lambda.ineqlin-lambda的线性不等式；
       lambda.eqlin-lambda的线性等式。
三． 举例
例1：求解线性规划问题：
              max f=2x1+5x2
              s.t
先将目标函数转化成最小值问题：min(-f)=- 2x1-5x2
程序：
f=[-2 -5];
A=[1 0;0 1;1 2];
b=[4;3;8];
[x,fval]=linprog(f,A,b)
f=fval*(-1)
结果：   x = 2 
3
                     fval = -19.0000
maxf =  19
例2：minf=5x1-x2+2x3+3x4-8x5
s.t  -2x1+x2-x3+x4-3x5≤6
    2x1+x2-x3+4x4+x5≤7
    0≤xj≤15  j=1,2,3,4,5
程序：
f=[5 -1 2 3 -8];
A=[-2 1 -1 1 -3;2 1 -1 4 1];
b=[6;7];
lb=[0 0 0 0 0];
ub=[15 15 15 15 15];
[x,fval]=linprog(f,A,b,[],[],lb,ub)
结果：x =
           0.0000
           0.0000
           8.0000
           0.0000
            15.0000
minf =
  -104
例3：求解线性规划问题：
                     minf=5x1+x2+2x3+3x4+x5
s.t  -2x1+x2-x3+x4-3x5≤1
                2x1+3x2-x3+2x4+x5≤-2
                    0≤xj≤1  j=1,2,3,4,5
程序：
       f=[5 1 2 3 1];
       A=[-2 1 -1 1 -3;2 3 -1 2 1];
       b=[1;-2];
       lb=[0 0 0 0 0];
       ub=[1 1 1 1 1];
       [x,fval,exitflag,output,lambda]=linprog(f,A,b,[],[],lb,ub)                           运行结果：        
       Exiting: One or more of the residuals, duality gap, or total relative error
         has grown 100000 times greater than its minimum value so far:
         the primal appears to be infeasible (and the dual unbounded).
         (The dual residual < TolFun=1.00e-008.)

x = 0.0000
                   0.0000
                   1.1987
                   0.0000
                    0.0000
fval =
                  2.3975
exitflag =
                  -1
output =
          iterations: 7
           cgiterations: 0
         algorithm: 'lipsol'
lambda =
                  ineqlin: [2x1 double]
                 eqlin: [0x1 double]
                 upper: [5x1 double]
                 lower: [5x1 double]
       显示的信息表明该问题无可行解。所给出的是对约束破坏最小的解。
       例4：求解实例1的生产计划问题
建立数学模型：
设x1、x2分别为生产甲、乙产品的件数。f为该厂所获总润。
       max f=70x1+120x2
       s.t  9x1+4x2≤3600
        4x1+5x2≤2000
        3x1+10x2≤3000
        x1,x2≥0
将其转换为标准形式：
min f=-70x1-120x2
       s.t  9x1+4x2≤3600
        4x1+5x2≤2000
        3x1+10x2≤3000
        x1,x2≥0

       程序：   f=[-70 -120];
                     A=[9 4 ;4 5;3 10 ];
                     b=[3600;2000;3000];
                     lb=[0 0];
                     ub=[];
                            [x,fval,exitflag]=linprog(f,A,b,[],[],lb,ub)
                            maxf=-fval
              结果：   x =
                           200.0000
                           240.0000
fval =
                                  -4.2800e+004
exitflag =
    1
maxf =
      4.2800e+004
      例5：求解实例2
       建立数学模型：
max f=0.15x1+0.1x2+0.08 x3+0.12 x4
       s.t  x1-x2- x3- x4≤0
        x2+ x3- x4≥0
              x1+x2+x3+ x4=1
              xj≥0  j=1,2,3,4   
将其转换为标准形式：
min z=-0.15x1-0.1x2-0.08 x3-0.12 x4
       s.t  x1-x2- x3- x4≤0
        -x2- x3+ x4≤0
              x1+x2+x3+ x4=1
              xj≥0  j=1,2,3,4
       程序：   f = [-0.15;-0.1;-0.08;-0.12];
A =  [1 -1 -1 -1
                          0 -1 -1 1];
b = [0; 0];
Aeq=[1 1 1 1];
beq=[1];
lb = zeros(4,1);
                    [x,fval,exitflag] = linprog(f,A,b,Aeq,beq,lb)
                     f=-fval
       结果：x =
                         0.5000
                         0.2500
                         0.0000
                         0.2500
fval =
                          -0.1300
exitflag =
                               1
f =
0.1300
       即4个项目的投资百分数分别为50%，25%，0,  25%时可使该公司获得最大的收益，其最大收益可到达13%。过程正常收敛。
       例6：求解实例3
       建立数学模型：
设ai j为由工厂i运到市场j的费用，xi j 是由工厂i运到市场j的箱数。bi是工厂i的产量，dj是市场j的需求量。

b= ( 60 40 50 )T   d= ( 20 35 33 34 )T

       s.t 

        x i j≥0
       程序：   A=[2 1 3 2;1 3 2 1;3 4 1 1];
                     f=A(:);
                     B=[ 1 0 0 1 0 0 1 0 0 1 0 0
                            0 1 0 0 1 0 0 1 0 0 1 0
                            0 0 1 0 0 1 0 0 1 0 0 1];
                     D=[1 1 1 0 0 0 0 0 0 0 0 0
                            0 0 0 1 1 1 0 0 0 0 0 0
                            0 0 0 0 0 0 1 1 1 0 0 0
                            0 0 0 0 0 0 0 0 0 1 1 1];
                     b=[60;40;50];
                     d=[20;35;33;34];
                     lb=zeros(12,1);
                     [x,fval,exitflag]=linprog(f,B,b,D,d,lb)
       结果：   x =
                         0.0000
                          20.0000
                         0.0000
                          35.0000
                         0.0000
                         0.0000
                         0.0000
                         0.0000
                          33.0000
                         0.0000
                          18.4682
                          15.5318
fval =
                            122.0000
exitflag =
     1
       即运输方案为：甲市场的货由B厂送20箱；乙市场的货由A厂送35箱；丙商场的货由C厂送33箱；丁市场的货由B厂送18箱，再由C厂送16箱。
最低总运费为：122元。

§2 非线性规划模型
一．非线性规划课题
实例1  表面积为36平方米的最大长方体体积。
建立数学模型：
设x、y、z分别为长方体的三个棱长，f为长方体体积。
max f = x y (36-2 x y)/2 (x+y)
实例2  投资决策问题
某公司准备用5000万元用于A、B两个项目的投资，设x1、x2分别表示配给项目A、B的投资。预计项目A、B的年收益分别为20%和16%。同时，投资后总的风险损失将随着总投资和单位投资的增加而增加，已知总的风险损失为2x12+x22+(x1+x2)2.问应如何分配资金，才能使期望的收益最大，同时使风险损失为最小。
建立数学模型：
       max f=20x1+16x2-λ[2x12+x22+(x1+x2)2]
       s.t  x1+x2≤5000
        x 1≥0,x2≥0
目标函数中的λ≥0是权重系数。
由以上实例去掉实际背景，其目标函数与约束条件至少有一处是非线性的，称其为非线性问题。
非线性规划问题可分为无约束问题和有约束问题。实例1为无约束问题，实例2为有约束问题。
二．无约束非线性规划问题：
求解无约束最优化问题的方法主要有两类：直接搜索法(Search method)和梯度法(Gradient method).
1．fminunc函数
调用格式： x=fminunc(fun,x0)
                     x=fminunc(fun,x0,options)
                     x=fminunc(fun,x0,options,P1,P2)
              [x,fval]=fminunc(…)
              [x,fval, exitflag]=fminunc(…)
              [x,fval, exitflag,output]=fminunc(…)
[x,fval, exitflag,output,grad]=fminunc(…)
[x,fval, exitflag,output,grad,hessian]=fminunc(…)
说明：fun为需最小化的目标函数，x0为给定的搜索的初始点。options指定优化参数。
返回的x为最优解向量；fval为x处的目标函数值；exitflag描述函数的输出条件；output返回优化信息；grad返回目标函数在x处的梯度。Hessian返回在x处目标函数的Hessian矩阵信息。
例1 ： 求
程序：编辑ff1.m文件
function f=ff1(x)
f=8*x(1)-4*x(2) +x(1)^2+3*x(2)^2;
通过绘图确定一个初始点：
[x,y]=meshgrid(-10:.5:10);
z= 8*x-4*y +x.^2+3*y.^2;
surf(x,y,z)
选初始点：x0=(0,0)
x0=[0,0];
[x,fval,exitflag]=fminunc(@ff1,x0)

结果：x =
            -4.0000    0.6667
fval =
             -17.3333
exitflag =
                 1
例2：
程序：编辑ff2.m文件：
function f=ff2(x)
f=4*x(1)^2+5*x(1)*x(2)+2*x(2)^2;
取初始点：x0=(1,1)
x0=[1,1];
[x,fval,exitflag]=fminunc(@ff2,x0)
       结果：    x =
                           1.0e-007 *
                          -0.1721    0.1896
fval =
                           2.7239e-016
exitflag =
                               1
       例3：将上例用提供的梯度g最小化函数进行优化计算。
修改M文件为：
function [f,g]=ff3(x)
f=4*x(1)^2+5*x(1)*x(2)+2*x(2)^2;
if nargut >1
     g(1)=8*x(1)+5*x(2);
     g(2)=5*x(1)+4*x(2);
end
通过下面将优化选项结构options.GradObj设置为’on’来得到梯度值。
       options=optimset(‘Gradobj’,’on’);
       x0=[1,1];
[x,fval,exitflag]=fminunc(@ff3,x0,options)
       结果：   x =
                           1.0e-015 *
                          -0.2220   -0.2220
fval =
                           5.4234e-031
exitflag =
                        1
2． minsearch函数
调用格式： x=fminsearch(fun,x0)
                     x=fminsearch(fun,x0,options)
                     x=fminsearch(fun,x0,options,P1,P2)
              [x,fval]=fminsearch(…)
              [x,fval, exitflag]=fminsearch(…)
              [x,fval, exitflag,output]=fminsearch(…)
[x,fval, exitflag,output,grad]=fminsearch(…)
[x,fval, exitflag,output,grad,hessian]=fminsearch(…)
说明：参数及返回变量同上一函数。对求解二次以上的问题，fminsearch函数比fminunc函数有效。    

3． 多元非线性最小二乘问题：
非线线性最小二乘问题的数学模型为：
其中L为常数。
调用格式：  x=lsqnonlin(fun,x0)
                     x=lsqnonlin(fun,x0,lb,ub)
x=lsqnonlin(fun,x0,options)
                     x=lsqnonlin(fun,x0,options,P1,P2)
              [x,resnorm]=lsqnonlin(…)
              [x,resnorm, residual,exitflag]=lsqnonlin(…)
              [x,resnorm, residual , exitflag,output]=lsqnonlin(…)
[x,resnorm, residual,exitflag, output,lambda]=lsqnonlin(…)
[x,resnorm, r esidual,exitflag, output,lambda,jacobian]=lsqnonlin(…)
说明：x返回解向量；resnorm返回x处残差的平方范数值：sum(fun(x).^2)；residual返回x处的残差值fun(x)；lambda返回包含x处拉格朗日乘子的结构参数；jacobian返回解x处的fun函数的雅可比矩阵。
lsqnonlin默认时选择大型优化算法。Lsqnonlin通过将options.LargeScale设置为’off’来作中型优化算法。其采用一维搜索法。
例4．求 minf=4(x2-x1)2+(x2-4)2 ，选择初始点x0(1,1)
程序：
f ='4*(x(2)-x(1))^2+(x(2)-4)^2'
[x,reshorm]=lsqnonlin(f,[1,1])
结果：   x =
                  3.9896    3.9912
reshorm =
                    5.0037e-009
例5：求 ，选择初始点x0(0.2,0.3)
求解：先编辑ff5.m文件：
              function f=ff5(x)
              k=1:10;
              f=2+2*k-exp(k*x(1))-exp(k*x(2));
然后作程序：x0=[0.2,0.3];
                     [x,resnorm]=lsqnonlin(@ff5,x0)
结果 ： x =
                  0.2578    0.2578
resnorm =
                    124.3622
二． 有约束非线性规划问题：
数学模型：  min F(x)
s.t  Gi (x) ≤0             i=1,…,m
                 Gj (x) =0         j=m+1,…,n
                 xl≤x≤xu
其中：F(x)为多元实值函数，G(x)为向量值函数，
在有约束非线性规划问题中，通常要将该问题转换为更简单的子问题，这些子问题可以求并作为迭代过程的基础。其基于K-T方程解的方法。它的K-T方程可表达为：

       方程第一行描述了目标函数和约束条件在解处梯度的取消。由于梯度取消，需要用拉格朗日乘子λi来平衡目标函数与约束梯度间大小的差异。
调用格式:    x=fmincon(f,x0,A,b)
                            x=fmincon(f,x0,A,b,Aeq,beq)
                            x=fmincon(f,x0,A,b,Aeq,beq,lb,ub)
                            x=fmincon(f,x0,A,b,Aeq,beq,lb,ub,nonlcon)
                            x=fmincon(f,x0,A,b,Aeq,beq,lb,ub,nonlcon,options)
              [x,fval]=fmincon(…)
              [x, fval, exitflag]=fmincon(…)
              [x, fval, exitflag, output]=fmincon(…)
              [x, fval, exitflag, output, lambda]=fmincon(…)
       说明：x=fmincon(f,x0,A,b)返回值x为最优解向量。其中：x0为初始点。A,b为不等式约束的系数矩阵和右端列向量。
       x=fmincon(f,x0,A,b,Aeq,beq) 作有等式约束的问题。若没有不等式约束，则令A=[ ]、b=[ ] 。
x=fmincon(f, x0,A,b,Aeq,beq,lb,ub, nonlcon ,options) 中lb ,ub为变量x的下界和上界；nonlcon=@fun,由M文件fun.m给定非线性不等式约束c (x) ≤0和等式约束g(x)=0；options为指定优化参数进行最小化。
例6：求解：min 100(x2-x12 )2+(1-x1)2
          s.t  x1≤2;
x2≤2
程序：首先建立ff6.m文件：
function f=ff6(x)
f=100*(x(2)-x(2)^2)^2+(1-x(1))^2;
然后在工作空间键入程序：
x0=[1.1,1.1];
A=[1 0;0 1];
b=[2;2];
[x,fval]=fmincon(@ff6,x0,A,b)
结果：  x =
                1.0000    1.0000
fval =
        3.1936e-011
    例7：求解：

        首先建立目标函数文件ff7.m文件：
            function f=ff7(x)
            f=-x(1)*x(2)*x(3)
        然后将约束条件改写成如下不等式：
                -x1-2x2-2x3≤0
                x1+2x2+2x3≤72
    在工作空间键入程序：
            A=[-1 -2 -2;1 2 2];
            b=[0;72];
            x0=[10;10;10];
            [x,fval]=fmincon(@ff71,x0,A,b)
    结果：  x =
                24.0000
                12.0000
                12.0000
fval =
                -3456
例8求解：minf=ex1(6x12+3x22+2x1x2+4x2+1)
              s.t   x1x2-x1-x2+1≤0
            -2x1x2-5≤0
程序：首先建立目标函数文件ff8.m文件：
        function  f=ff8(x)
              f=exp(x(1))*(6*x(1)^2+3*x(2)^2+2*x(1)*x(2)+4*x(2)+1);
       再建立非线性的约束条件文件：ff8g.m
              function [c,g]=ff8g(x)
              c(1)=x(1)*x(2)-x(1)-x(2)+1；
              c(2)=-2*x(1)*x(2)-5；
g=[];
然后在工作空间键入程序：
        x0=[1,1];
              nonlcon=@ff8g
[x, fval] =fmincon(@ff8,x0,[],[],[],[],[],[], nonlcon)
结果：  x =
            -2.5000    1.0000
fval =
            3.3244
exitflag =
            1
    当有等式约束时，要放在矩阵g的位置，如上例中加等式约束：
        x(1)+2*x(1)=0
程序：首先建立 fun1.m文件：
        function[c,g]=ff8g1(x)
       c(1)=x(1)*x(2)-x(1)-x(2)+1；
       c(2)=-2*x(1)*x(2)-5；
       g(1)=x(1)+2*x(2);
然后在工作空间键入程序：
        x0=[-1,1];
       nonlcon=@ff8g1;
        [x, fval,exitflag] =fmincon(@ff8,x0,[],[],[],[],[],[], nonlcon)
结果：  x =
                -2.2361    1.1180
fval =
                3.6576
exitflag =