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