建模训练案例-2：简单线性规划问题

作者：祁彬彬、马良
本篇问题引自《数学模型(第五版)》
序

  本篇文章主要讲解很多教材在开篇时，都会介绍的一个线性规划求解的简单实例——奶制品加工生产计划获利的代码，但它出现在我们关于建模文章系列的开篇位置，其承担的任务有所不同，因为写出代码方案不再是目标，而是手段。我们关注的是，通过多种求解代码方案，几个非常关键的问题应当在这个系列刚开始的时候就得到部分回答：

• Lingo作为早期数学建模中的应知必会软件，其模型搭建与求解机制是怎样的？
• 为什么早期版本的MATLAB在优化问题求解时无法代替Lingo？
• 为什么新版本MATLAB又具备了走向建模舞台中央，彻底取代Lingo的资格和条件？
  

  这几个问题，尤其第3个，兼具头条系的帅气，和UC震惊部的沧桑，究竟是博眼球？还是确有资格？可能需要一个比较漫长的论证和测试，并伴随大量问题实例以及对应的代码求解方案。

那么从这篇文章开始，我们可能就要开始这个探讨过程，这篇文章主要有如下内容：

• 用Lingo软件求解该优化模型；
• 分析总结Lingo代码中，有关优化模型求解时对应的基本要素与机制；
• 利用Yalmip求解同一问题的模型构造与求解方法；
• 利用MATLAB求解同一问题，其对应Lingo不同模块位置的命令函数以及相应的代码编写方法；
• 设置MATLAB调用Gurobi求解问题的代码方案与步骤。
  

通过对这个示例的剖析，达到如下目标：

• 解释Lingo模型搭建与求解机制以及新版MATLAB中，对应的等效代码流程；
• 通过比较几种方式求解同一规划问题的代码，感受MATLAB基于问题和基于求解器两类求解代码的编写方式差异；
• 了解MATLAB调用Gurobi求解同一模型的方法，以及如何在代码中指定调用MATLAB官方命令或者Gurobi同名函数求解模型；
• 利用Yalmip建模并求解(调用Gurobi)的方式；
  

题目  一奶制品加工厂用牛奶生产​ 两种奶制品，1桶牛奶可在甲车间用12h加工成3kg的​ ，或者在乙车间用8h加工成4kg的 ​，根据市场需求，生产出的​ 全部都能售出，​每千克 获利24元，​每千克 获利16元，现在加工厂每天能得到50桶牛奶的供应，每天正式工人总的劳动时间为480h，且甲车间设备每天至多能加工100kg的​ ，乙车间的设备加工能力可以认为没有上限限制(加工能力足够大)，试为该厂指定一个生产计划，使得每天的获利最大。

分析

获利最大的决策是生产计划，每天多少桶牛奶用来各自生产 ​，决策受3个条件限制：原料供应、劳动时间、甲车间加工能力。

优化模型要素

• 决策变量：设每天用 ​桶牛奶生产 ​、用​ 桶牛奶生产 ​;
• 目标函数：设每天获利 ​元。​ 桶牛奶生产​ 的 ​，获利为​ ，​ 桶牛奶生产​ 的 ​，获利为​ ，所以​
• 约束条件：包括原料供应、劳动时间、设备能力、非负约束。
  

约束条件明细：

• 原料供应：生产总量所需牛奶上限需满足​   ：
• 劳动时间：加工总时间不超过工人总劳动时间  ​：
• 设备能力：加工的产量上限   ​：
• 非负约束  
  

数学模型

从优化模型的几个要素，以及数学模型表达式，很容易看出牛奶加工问题属于简单线性规划模型求解，这种优化模型任何市面所见的优化软件都不会有求解难度，且每一条约束和目标函数都容易理解，后续分析的重心放在代码的流程结构上。

代码Lingo与MATLAB求解器求解模型方式分析比较Lingo求解优化问题的机制

一般优化问题从模型搭建到求解，要经历下面的过程。

  需要强调的是，建模语言进入求解器之前，总归要被转换为矩阵形式输入的，但Lingo却中间两个环节做了整合化的处理，不妨看看Lingo是如何求解奶制品加工厂这个简单LP例子的：

1. max = 72*x1+64*x2;
   
2. x1+x2 <= 50;
   
3. 12*x1+8*x2 <= 480;
   
4. 3*x1<=100;
   
5. <span style="color: rgb(68, 68, 68); font-family: Tahoma, &quot;Microsoft Yahei&quot;, Simsun; font-size: 14px; font-style: normal; font-variant-ligatures: normal; font-variant-caps: normal; font-weight: 400; letter-spacing: normal; orphans: 2; text-align: start; text-indent: 0px; text-transform: none; white-space: normal; widows: 2; word-spacing: 0px; -webkit-text-stroke-width: 0px; background-color: rgb(255, 255, 255); text-decoration-thickness: initial; text-decoration-style: initial; text-decoration-color: initial; display: inline !important; float: none;"></span>

复制代码

  和数学模型比较，会发现Lingo代码几乎一模一样，不过，这不是计算机能够读懂的模型形式，但它符合人的数学思维——用户只需要把优化模型“翻译”到Lingo里就可以了，Lingo内部会在求解前把“用户优化模型”转换为“计算机优化模型”，最后为用户返回如下计算结果：

1. Global optimal solution found.
   
2.   Objective value:                              3360.000
   
3.   Infeasibilities:                                  0.000000
   
4.   Total solver iterations:                             2
   
5.   Elapsed runtime seconds:                        0.07
   
6. 
   
7.   Model Class:                                        LP
   
8. 
   
9.   Total variables:                          2
   
10.   Nonlinear variables:                  0
    
11.   Integer variables:                      0
    
12. 
    
13.   Total constraints:                       4
    
14.   Nonlinear constraints:               0
    
15. 
    
16.   Total nonzeros:                         7
    
17.   Nonlinear nonzeros:                 0
    
18. 
    
19. 
    
20. 
    
21.                                 Variable         Value        Reduced Cost
    
22.                                       X1          20.00000            0.000000
    
23.                                       X2          30.00000            0.000000
    
24. 
    
25.                                      Row    Slack or Surplus      Dual Price
    
26.                                        1            3360.000            1.000000
    
27.                                        2            0.000000            48.00000
    
28.                                        3            0.000000            2.000000
    
29.                                        4            40.00000            0.000000

复制代码

自第21行向下即为Lingo的运行结果：当前条件下，购买20桶牛奶加工​、30桶牛奶加工​，这个最优生产计划可以得到最大利润3360元。

从上述问题的求解可以看出，Lingo在一定程度上，封装了优化模型中，与计算机对接的模型转换构造部分。

牺牲了一点点转换效率，却大幅提高了用户模型搭建，以及学习Lingo表达优化模型的效率，这就是Lingo长期以来在高校数学建模比赛中，深受学生和老师欢迎，且长盛不衰的原因之一。

进一步地，Lingo还可以把上述模型更加参数化和标准化，以变量定义“sets”、数据输入“data”、数据计算“cals”、处置设置“inti”和约束“[Obj & Con]”这5个部分，重写前述程序：

1. sets:
   
2. milk/1..2/:X, bulk_to_kg, price,work_time;
   
3. endsets
   
4. 
   
5. data:
   
6. bulk_to_kg = 3,4;
   
7. price = 24, 16;
   
8. work_time = 12, 8;
   
9. total_milk = 50;
   
10. upper_All = 100;
    
11. work_time_all = 480;
    
12. enddata
    
13. 
    
14. calc:
    
15. endcalc
    
16. 
    
17. init:
    
18. endinit
    
19. 
    
20. [objective]max = @sum(milk(i):X(i)*bulk_to_kg(i)*price(i));
    
21. [raw]@sum(milk(i):X(i)) <= total_milk;
    
22. [work_time_cons]@sum(milk(i):work_time(i)*X(i)) <= work_time_all;
    
23. [instrument]bulk_to_kg(1)*X(1) <= upper_All;

复制代码

运行结果：

1. Global optimal solution found.
   
2.   Objective value:                              3360.000
   
3. ...
   
4.                                 Variable                     Value            Reduced Cost
   
5.                               TOTAL_MILK             50.00000             0.000000
   
6.                                UPPER_ALL            100.0000               0.000000
   
7.                            WORK_TIME_ALL       480.0000               0.000000
   
8.                                    X( 1)                     20.00000             0.000000
   
9.                                    X( 2)                     30.00000             0.000000
   
10.                           BULK_TO_KG( 1)             3.000000           0.000000
    
11.                           BULK_TO_KG( 2)             4.000000           0.000000
    
12.                                PRICE( 1)                 24.00000              0.000000
    
13.                                PRICE( 2)                 16.00000              0.000000
    
14.                            WORK_TIME( 1)          12.00000              0.000000
    
15.                            WORK_TIME( 2)            8.000000            0.000000
    
16. 
    
17.                                      Row               Slack or Surplus           Dual Price
    
18.                                OBJECTIVE               3360.000                   1.000000
    
19.                                      RAW                        0.000000            48.00000
    
20.                           WORK_TIME_CONS            0.000000              2.000000
    
21.                               INSTRUMENT                40.00000               0.000000

复制代码

这个例子没有用到calc和inti这两个模块的设置，计算结果与之前相同。
MATLAB基于求解器建模的机制  看完Lingo求解问题的模式，接下来再看MATLAB的做法，它是基于矩阵思维的程序编写语言，这个特色使得它几十年以来长期保持着工程计算软件领域的先锋地位，因此在MATLAB介入优化计算领域伊始，优化工具箱命令譬如：linprog、fminunc、fmincon等，都秉承着矩阵化的函数设计思路，为了说明这种风格的代码形式，不妨以linprog再度求解上述奶制品生产计划LP问题。

1. >> [x,fval] = linprog(-[72 64],[1 1;12 8;3 0],[50 480 100]',[],[],[0 0],[300 300])
   
2. Optimal solution found.
   
3. x =
   
4.    20.0000
   
5.    30.0000
   
6. fval =
   
7.        -3360

复制代码

  代码仅有一行，甚至比Lingo更简洁，结果也和Lingo相同，不过对于初学者，代码还是有必要解释一下，linprog的调用格式如下：

1. [x, fvl] = linprog(f,A,b,Aeq,beq,lb,ub)

复制代码

  MATLAB中，lingprog有多种调用方式，上述是和本例调用完全对应的一种，如果把原来的代码稍微“解压缩”，会更容易看出基于求解器的模型编写特点：

1. f = -[72, 64];                % 矩阵表示的目标函数
   
2. 
   
3. A = [1 1; 12 8; 3 0];     % 矩阵表示的不等式约束的左侧表达式
   
4. b = [50 480 100]';       % 矩阵表示的不等式约束的右侧系数向量
   
5. 
   
6. Aeq = [];                     % 矩阵表示的等式约束左侧表达式
   
7. beq = [];                      % 矩阵表示的等式约束右侧系数向量
   
8. 
   
9. lb = [0 0];                    % 设计变量下界
   
10. ub = [300 300];           % 设计变量上界
    
11. 
    
12. [x, fvl] = linprog(f,A,b,Aeq,beq,lb,ub)  % 调用linprog求解

复制代码

上述代码的注释部分表明了各个参数的意义，观察MATLAB基于求解器的代码，会发现这样几个明显的特点：

• 无论目标函数或者约束条件，代码在定义时，一律基于矩阵形式构造，例如第1行代码“f = -[72 64]”代表目标函数为：​，而模型中的三组形式类似的不等式约束条件则分别用第2,3参数​，写成如下形式，用于代表“​”：
  

                                                         

• 即使本例中没有等式约束，要用空矩阵“[]”占位；
• 整个代码的模型搭建过程中，优化模型的设计变量似乎是“虚”的，因为既没有定义设计变量​，形式上当然也没出现在矩阵化模型当中；
• 模型只能求解目标函数的最小值，因此奶制品生产计划的利益最大化，需要在目标函数中利用负号手动处理；
• 约束条件只接受以小于等于形式，因此如果有带大于的条件，也需要用负号手动处理；
• 调用linprog的过程是函数体形式的，也就是函数名被明确指定为linprog。
  
  

比较Lingo和MATLAB基于求解器式的代码求解方式，表面上似乎看不出什么优劣高下，貌似无非是两种不同的软件根据各自风格特点，形成了迥异的模型搭建方式而已，何况结果都是相同的。

实际完全不是这么回事儿。

MATLAB基于求解器的优化模型构造，从数学模型到矩阵形式，需用户去手动解析，这样的话，简单问题虽然要花一些功夫，但仍然可以解决，比如本例奶制品加工计划模型即循此例；不过，当优化模型约束条件趋于复杂化，数学形式到矩阵的模型抽象和解析过程，可能将花费用户大量的精力，而且许多错误比较隐蔽不易排查——试想一下，即使一个上千条约束规模的问题，转换成矩阵形式时，由于不同类型的约束条件具有不同的维度，因此就要对原约束条件被动扩维，然而扩维之后在第579个条件的第849列数据处，不小心发生了一个输入错误...

但产生这类麻烦的可能性，Lingo中从一开始就被降低到可被用户接受的程度，因为Lingo约束条件的写流程倾向于“碎片化”逐条输入，没有被动扩维，条件中变量意义明晰，约等于直接输入，无需在诸如：“​”这样条件里的​和​之间，扩充很多个0，而反观MATLAB，即使拥有稀疏矩阵求解函数sparse，但这类约束条件的写法仍然非常麻烦，加上自r2014b版本才正式拥有求解线性整数规划的求解能力(intlinprog)、且全局规划求解工具箱一直颇有些捉急的计算效率，一系列历史原因，导致MATLAB优化求解能力并没有在数学建模各类竞赛中，获得不俗的口碑。只是因为强悍的数据处理以及绘图能力，才让MATLAB没有淡出数学建模的舞台，但却始终不温不火，扮演着“第二火力输出点”的角色。MATLAB+Yalmip求解


  针对MATLAB一直以来在复杂模型表述方面的弱点与短板，瑞士林雪平大学的Johan Löfberg开发了第三方的Yalmip工具箱，可以通过问题式建模形式实现优化模型的描述以及求解，例如本例如果用Yalmip搭建并求解优化模型问题，其完整代码如下：

1. %% 1.sets
   
2. n = 2;
   
3. X = sdpvar(1,n);
   
4. %% 2.data
   
5. bulk_to_kg = [3,4];
   
6. price = [24, 16];
   
7. work_time = [12, 8];
   
8. total_milk = 50;
   
9. upper_A1 = 100;
   
10. work_time_all = 480;
    
11. 
    
12. %% 5.object_and_constraints
    
13. % f = X*bulk_to_kg'*price'; 向量操作，但对复杂分支结构未必通用
    
14. % c1 = sum(X)<=total_milk;
    
15. 
    
16. f = 0;
    
17. for i = 1 : n
    
18.     f = f + X(i)*bulk_to_kg(i)*price(i);
    
19. end
    
20. [c1, c2] = deal(0);
    
21. for i = 1 : n
    
22.     c1 = c1 + X(i);
    
23.     c2 = c2 + work_time(i)*X(i);
    
24. end
    
25. c3 = 3 * X(1);
    
26. C = [c1<=total_milk, c2<=work_time_all,c3<=upper_A1];
    
27. 
    
28. % solve
    
29. options = sdpsettings('solver','gurobi','verbose', 0);
    
30. diagnostics  = optimize(C,-f,options);
    
31. x = value(X)
    
32. o = -value(f)

复制代码

上述代码标出了与Lingo对应的sets、data和obj & constraints对等部分，求解时，用sdpsettings指定gruobi为求解器，有两点值得注意：

• Yalmip是以MATLAB为平台，再调用第三方求解器Gurobi的；
• Yalmip以MATLAB为平台，需要遵从MATLAB的基本语法流程，事实上它自身就是在MATLAB早期的对象化编程体系下编写出来的工具箱，因此在Yalmip中编写优化模型，既可以借助MATLAB自身编程体系中的优势，也能够享受Yalmip问题式建模的红利。
  
  

注：上述代码第`12~14`行的注释部分，是下方循环体内目标函数`f`和约束`c1`的向量化写法，二者在约束的表示效果方面是完全相同的，这里采用循环一方面是为对应Lingo代码，便于理解；此外也说明MATLAB在表述同一种问题时，所能采取编程手段的多样性；最后，如果初学MATLAB，能在建模的代码中随手用出矢量化的代码，可能需要一段时间训练，但在后续代码中，我们将在优化实例中越来越频繁地采用矢量化方式的代码，也期望经过我们文章中实例代码的洗礼，在今后建模比赛实际论文中，看到越来越多的矢量化代码应用。

因此总结下来，Yalmip工具箱的好处显而易见，主要表现在：

• 更多的求解器接口：它提供市面上绝大多数优化求解器的接口，调用方法标准化且代码非常简洁，在上述示例求解奶制品生产计划问题时，以MATLAB为平台，通过Yalmip调用第三方求解器Gurobi就是一个典型案例；
• 更兼容的语法：以MATLAB为平台，却兼具通过基于问题模优化模型描述及MATLAB自身程序编写语法支持两大优点，这让非常复杂的约束条件编写时，代码代价更低，不但能以类似Lingo的问题式流程描述模型，还拥有Lingo不具备的MATLAB语法背景支持；
• 更丰富的求解器支持：Yalmip提供一些MATLAB至今仍然不具备的求解器运算功能，例如分支定界算法求解非线性混合整数规划问题，这一点在后续本系列文章中会陆续讨论，因此，以Yalmip作为中介环境，MATLAB中编写优化模型代码并计算，堪称是颠覆式的全新优化选择。
  
  

MATLAB问题式建模

  MATLAB在r2017b版本推出“问题式”优化模型描述方式，这是MATLAB在全新对象化编程体系下，为优化工具箱订制的模型描述架构，基于独立优化对象基类optim，最大程度复现问题式建模从设计变量(optimvar)到问题与约束条件的定义和描述(optimproblem)。且求解器可通过MATLAB自行指定，这实际上把linprog、intlinprog、fmincon、quadprog、lsqnonlin等内置优化求解器纳入一个完整而标准化的体系中，使得用户可以集中精力去描述优化模型，简化了求解流程。

仍以本文的奶制品生产计划模型为例，对这类简单的线性规划问题，前面已经说过MATLAB调用的函数是linprog，但以何种形式调用linprog函数，是个需要重点关注的问题。自r2017b版本起，MATLAB增添了“基于问题的模型构造方式”(在r2019b再次对问题式建模的代码流程做了调整，增加了几个类似补丁的函数，这种差异在后续文章中再结合实例讲解)，这是更适合于描述复杂优化模型的代码书写方式，下面扼要讲解新版本的模型求解代码。

模型代码与运行结果首先，按照数学模型，可以写出如下代码：

1. x = optimvar('x',2, 'LowerBound',0);                               % 设计变量
   
2. prob = optimproblem('ObjectiveSense','max','Objective',[72 64]*x); % 决策目标
   
3. prob.Constraints.con = [1 1;12 8;3 0]*x<=[50;480;100];             % 约束条件
   
4. [sol,fvl] = solve(prob)
   
5. sol.x

复制代码

运行结果：

1. Solving problem using linprog.
   
2. Optimal solution found.
   
3. sol =
   
4.   struct with fields:
   
5.     x: [2×1 double]
   
6. fvl =
   
7.         3360
   
8. >> sol.x
   
9. ans =
   
10.    20.0000
    
11.    30.0000

复制代码

上述结果仍然是：加工厂每天用20桶牛奶生产​、用30桶牛奶生产​，可以得到最大收益3360元，结果和Lingo、或者Yalmip调用Gurobi完全相同，但这里我们还需要额外讨论一下MATLAB问题式建模的代码流程。
MATLAB问题式建模代码剖析

注意MATLAB基于问题建模的代码是如何表达优化设计三要素的：

• 优化变量：函数optimvar创建​实型优化变量，指定设计变量下界均为0，即：牛奶桶数大于零。如果查看工作空间，会发现，生成的设计变量x，其类型的名称为optim.problemdef.OptimizationVariable；
• 优化模型：函数optimproblem创建基于问题描述方式的优化模型，如果查看工作空间，会发现所产生的的问题模型，其类型为optim.problemdef.OptimizationProblem；
  
  

MATLAB初学者也许还没有关注变量类型的习惯，但对于有一定MATLAB使用经验的用户，会困惑于：数据类型为什么看起来这么复杂？实际上，问题式优化模型描述方式，主要基于面向对象的编程逻辑：在类型名称里的点调用，不是结构数组，而是属于“类”的定义描述范畴。例如：optim.problemdef.OptimizationProblem，optim是优化模型描述的父类，problemdef.OptimzationProblem是继承父类的子类，所生成的具体模型prob，则是基于优化问题定义类的一种实例对象，它通过参数，具体定义了目标函数求取最大目标值，并给出确定的目标函数[72,64]*x，这就把原先抽象的类，借助模型数据把行为给具体化了

注：以上知识属于拓展，如果不理解对象化的思想也没关系，它的掌握与否，不会影响建模代码的编写，我们只需要知道，无论设计变量的类型如何奇怪，但它既支持矩阵表述，也可以作为普通形式的变量，在整个问题描述和求解过程中，被当做普通变量参与运算就可以了。

  定义了目标函数、设计变量，接着是约束条件，它同样是基于实例化问题prob的子类，注意前述模型代码第3行，约束条件的名称为prob.Constraints.con，最后的con是这组约束条件的名称，由用户自己定义。在牛奶加工的问题中，约束条件只有一组，可以通过矩阵书写形式，很方便地合并为1组约束条件——这也是用MATLAB编写优化模型的优势之一。但如果有很多组约束条件，无法通过矩阵或其他形式，统一放在一个con约束条件里，就可以在prob.Constraints下继续定义更多的con1、con2、...，例如：

1. prob.Constraints.con1 = ...
   
2. prob.Constraints.con2 = ...
   
3. ...

复制代码

MATLAB问题式模型的显示  当模型变得比较复杂，即使很仔细检查也不一定能完全写对，为提高模型的检查准确率，自r2019b新增show函数，可以把问题式模型还原成符合数学逻辑的形式：

1. >> show(prob)
   
2.   OptimizationProblem :
   
3.         Solve for:
   
4.        x
   
5.         maximize :
   
6.        72*x(1) + 64*x(2)
   
7.         subject to con:
   
8.        x(1) + x(2) <= 50
   
9.        12*x(1) + 8*x(2) <= 480
   
10.        3*x(1) <= 100
    
11.         variable bounds:
    
12.        0 <= x(1)
    
13.        0 <= x(2)

复制代码

这样对照原数学模型，就很容易检查出，模型到底写得对不对了。

为说明MATLAB问题式建模的模型描述机制，我们不妨也重写它的代码，使流程和Yalmip一样，也和Lingo对应起来，找到与sets、data和obj & constraints等对应的MATLAB代码：

1. %% 1.sets
   
2. n = 2;
   
3. X =  optimvar('X',n);
   
4. 
   
5. %% 2.data
   
6. bulk_to_kg = [3,4];
   
7. price = [24, 16];
   
8. work_time = [12, 8];
   
9. total_milk = 50;
   
10. upper_A1 = 100;
    
11. work_time_all = 480;
    
12. 
    
13. %% 5.object_and_constraints
    
14. % f = X*bulk_to_kg'*price'; 向量操作，但对复杂分支结构未必通用
    
15. % c1 = sum(X)<=total_milk;
    
16. 
    
17. prob = optimproblem('ObjectiveSense','maximize');
    
18. f = 0;
    
19. for i = 1 : n
    
20.     f = f + X(i)*bulk_to_kg(i)*price(i);
    
21. end
    
22. [c1, c2] = deal(0);
    
23. for i = 1 : n
    
24.     c1 = c1 + X(i);
    
25.     c2 = c2 + work_time(i)*X(i);
    
26. end
    
27. c3 = 3 * X(1);
    
28. 
    
29. prob.Objective = f;
    
30. 
    
31. prob.Constraints.raw_cons = c1<=total_milk;
    
32. prob.Constraints.worktime_cons = c2<=work_time_all;
    
33. prob.Constraints.instrument = c3<=upper_A1;
    
34. 
    
35. %% solve
    
36. [sol, fvl] = solve(prob);

复制代码

  观察上述代码，MATLAB新的问题式模型描述方式已经与Lingo在流程结构上具有很高的相似度，除了之前提到的问题式建模本身的优点，MATLAB官方问题式模型描述不但支持求解器模式的矩阵化约束条件表述，还能够实现一些求解器模式难以实现的细节，例如通过ObjectiveSense定义目标函数的最大化，或者直接采用“>=”定义大于等于约束，结合MATLAB的丰富函数库，能实现非常复杂的目标函数或约束条件的书写与订制，大幅度降低了用户的代码负担；最后，虽然截止目前，MATLAB官方的模型在指定求解器方面还比较封闭，但也并不是绝对与外部求解器隔绝，除了前面已经提到的，采用Yalmip构造模型并指定外部求解求解之外，当外部求解器与MATLAB具有良好的契合度，也可以通过一些比较另类的方式调用外部求解器，例如Gurobi。

1. >> x = optimvar('x',2, 'LowerBound',0);
   
2. prob = optimproblem('ObjectiveSense','max','Objective',[72 64]*x);
   
3. prob.Constraints.con = [1 1;12 8;3 0]*x<=[50;480;100];
   
4. addpath('C:\gurobi911\win64\examples\matlab')
   
5. [sol,fvl] = solve(prob)
   
6. 
   
7. Solving problem using linprog.
   
8. Academic license - for non-commercial use only - expires 2021-××-××
   
9. Gurobi Optimizer version 9.1.1 build v9.1.1rc0 (win64)
   
10. Thread count: × physical cores, × logical processors, using up to × threads
    
11. Optimize a model with 3 rows, 2 columns and 5 nonzeros
    
12. Model fingerprint: ×××××××××
    
13. Coefficient statistics:
    
14.   Matrix range     [1e+00, 1e+01]
    
15.   Objective range  [6e+01, 7e+01]
    
16.   Bounds range     [0e+00, 0e+00]
    
17.   RHS range        [5e+01, 5e+02]
    
18. Presolve removed 1 rows and 0 columns
    
19. Presolve time: 0.06s
    
20. Presolved: 2 rows, 2 columns, 4 nonzeros
    
21. Iteration    Objective       Primal Inf.    Dual Inf.      Time
    
22.        0   -3.6000000e+03   2.416667e+01   0.000000e+00      0s
    
23.        2   -3.3600000e+03   0.000000e+00   0.000000e+00      0s
    
24. Solved in 2 iterations and 0.12 seconds
    
25. Optimal objective -3.360000000e+03
    
26. sol =
    
27.   struct with fields:
    
28.     x: [2×1 double]
    
29. fvl =
    
30.         3360

复制代码

  当按照前一篇文章或官网的Lic注册要求的步骤，成功安装并设置，例如本机将Gurobi9.11以默认路径安装并设置好，利用addpath/rmpath这一对MATLAB函数手动添加/移除Gurobi内的同名MATLAB函数搜索路径，就可以实现直接调用/不调用Gurobi求解MATLAB模型问题了。

注：求解结果起始部分中，给出了求解器的基本信息，已经指明调用的是处于`Gurobi`安装文件夹下的`linprog`命令，而不是MATLAB官方函数`linprog`。

延伸讨论

根据前述模型以及求解代码，还可以继续分析如下几个问题

• 若35元可买到1桶牛奶，是否做这项投资？若投资，每天最多买多少桶牛奶？
• 若可以聘用临时工人以增加劳动时间，付给临时工人的工资最多是每小时多少元？
• 由于市场需求变化，每千克​获利增加到30元，是否应该改变生产计划？
  

延伸讨论1

原模型求解的最优结果为甲乙车间分别20和30桶牛奶，如果代入约束条件，发现牛奶原料被全部用完(​)，工人工时也全部用完(​)，但甲车间的生产力却并没有完全利用，这说明原料和工时属于“紧”约束，而甲车间的生产力却不是“紧”约束。

可以推知：如果紧约束被放松，效益必然还要增加，这样，效益就具有了潜在价值，这在经济上被称为“影子价格”，一桶牛奶的影子价格是48元，可是原料价格却只有35块，所以目前的情况下，如果加大牛奶原料的供应量，利润必然上升(13元/桶)，这自然是值得的投资。

如以桶为单位，增加牛奶供货量，在增至60桶前，利润都增加，意味着牛奶桶总数​时，影子价格都存在。影子价格在MATLAB中也可以直接求解，linprog命令中的第5个返回参数lambda返回解向量​的非零lagrange乘子系数，这里面包含的就是对偶价格：

1. x = optimvar('x',2, 'LowerBound',0); % 构造设计变量
   
2. prob = optimproblem('ObjectiveSense','max','Objective',[72 64]*x); % 构造决策目标
   
3. prob.Constraints.con = [1 1;12 8;3 0]*x<=[50;480;100];  % 构造约束条件
   
4. [sol,fvl,~,~,lam1] = solve(prob);
   
5. % -----------------------------------------------
   
6. >> lam1.Constraints.con
   
7. ans =
   
8.     48
   
9.      2
   
10.      0

复制代码

  第1个约束即为牛奶数量上限，可以看出其对偶价格就是之前分析的48，同理，第2约束为工人劳动时间的对偶价格，也同样存在一定效益提高上限，但提高空间不大。
延伸讨论2

与前一讨论类似，工人工时总量也是紧约束，如果增加一个工人工时，总利润会增加2元，所以如果聘用临时工，则付给工人的薪金不能超过2元/小时，且增加工人工时的上限是53.333h。

需要指出的是，MATLAB可以求解设计变量的影子价格，但它却不能像Lingo一样做敏感性分析，这项分析在某些场合下可能是会被用到的，因此我们可以通过软件间的协同，用MATLAB建模、借助Gurobi生成mps文件，由Lingo读取该mps文件，实现三软件的接力协同。

• 步骤1：我们自行编写了将MATLAB问题式优化模型转换为Gurobi模型的M程序，这实际上是按照规定的字段从MATLAB问题模型中读取相关的结构体域数据；
  
  

1. function model = prob2gurobi(prob)
   
2. model_s = prob2struct(prob);
   
3. model.A = [model_s.Aeq;model_s.Aineq];
   
4. model.obj = model_s.f;
   
5. model.rhs = [model_s.beq;model_s.bineq];
   
6. model.sense = repelem('=<',[numel(model_s.beq),numel(model_s.bineq)]);
   
7. model.vtype = repmat('C',1,length(model_s.lb));
   
8. model.vtype(model_s.intcon) = 'I';
   
9. model.modelsense = 'min';
   
10. model.lb    = model_s.lb;
    
11. model.ub    = model_s.ub;
    
12. end

复制代码

• 步骤2：在MATLAB中调用gurobi命令gurobi_write，把前述步骤得到的模型转换为mps文件；
  
  

1. model = prob2gurobi(prob);
   
2. gurobi_write(model, 'ExMilk.mps');

复制代码

• 步骤3：在Lingo中打开刚才的mps文件实现敏感性分析，这里需要在Lingo的Solver→Options→Dual Computations中，把默认的"Prices"改成“Prices & Ranges”,这样就可以在优化计算模型之后，按下Ctrl+R或Solver→Range查看敏感性分析报告了。运行结果：
  
  

1. Ranges in which the basis is unchanged:
   
2. 
   
3.                                        Objective Coefficient Ranges:
   
4. 
   
5.                                         Current        Allowable        Allowable
   
6.                       Variable      Coefficient         Increase         Decrease
   
7.                             C0         72.00000         24.00000         8.000000
   
8.                             C1         64.00000         8.000000         16.00000
   
9. 
   
10.                                            Righthand Side Ranges:
    
11. 
    
12.                                         Current        Allowable        Allowable
    
13.                            Row              RHS         Increase         Decrease
    
14.                             R0         50.00000         10.00000         6.666667
    
15.                             R1         480.0000         53.33333         80.00000
    
16.                             R2         100.0000         INFINITY         40.00000

复制代码

在上述结果中，可以看出设计变量数据的变化对最终利润，即：优化结果的影响，例如上述第15行说明，当当前正式工人的总劳动时间为480h，而这个值上浮上限(Allowable Increase)为53.3333h，也就是劳动时间按付给工人不超过￥2/h，继续增加至480+53.3333h这个阶段，是确定有利可图的。

采用上述方法得到敏感性分析报告，也应当注意如下几点：

• 三软件协同虽然看起来很酷，但实际上多数情况下没有必要，这样做是一些特定场合，可能要分析模型设计变量或约束条件的敏感性特征，而在MATLAB中求解，无法看到这种具体报告；
• 以上方法还可能用在另一个场景中，即：不熟悉Lingo，在Lingo中编写模型比较麻烦的情况；
• 上述方法并非通用，它受到了几个限制：首先，由于Gurobi不支持求解非线性整数规划模型，此类模型无法转换；其次，它仅支持对于普通linprog命令能够求解的LP模型。
  
  

延伸讨论3
如果每千克​获利增加到30元，代入原模型，把目标函数换成​即可，计算出的设计结果仍然是甲车间20桶，乙车间30桶，只是总利润上升到3720元，因此生产计划不用发生改变。这个结果也可以通过上述敏感性分析得到同样的结论，从上述敏感性分析报告第7行可以看出，第1个设计变量，也就是生产​的牛奶桶数，当它在​，也就是​范围内，模型基本结构，也就是生产计划无需改变就能获利。当​的获利从24元上升到30元，且​每桶牛奶生产3kg，单位获利数据是​，这个数据处于​允许范围之内，因此整个模型不需要重新计算，是必然获利的。总结

通过一道简单的LP问题，采用多种代码方式对其求解，主要目的是理解几个很关键的概念：

• MATLAB、Lingo等软件在求解优化问题时，模型搭建的主要机制与流程是可以被标准化的，这个标准化流程在不同优化求解器中是基本相同的；
• Lingo封装了以矩阵形式向求解器内部输入模型的过程，这种基于问题的建模方式是推荐方式；
• MATLAB早期版本无法实现基于问题建模的求解模式，而是被称为“基于求解器”式的建模形式，它需要手动把模型的目标函数、约束等，解析为矩阵模式才能求解，在复杂模型中，用户可能会为模型建立花费大量时间精力，并不是一种合理的模型搭建方式；
• Yalmip作为MATLAB环境下的第三方工具箱，基本弥补了MATLAB问题式建模的短板，甚至还可以再MATLAB环境下直接调用第三方求解器，或者还自带了基于分支定界算法的，可以进行混合整数非线性规划求解的求解函数，如果使用早期版本，它是利用MATLAB作为优化计算工具的强有力补充；
• MATLAB在r2017b推出了自己的问题式建模官方流程架构，这是基于全新的面向对象结构，在一个optim的基类下，可以流畅地使用MATLAB语言实现优化模型的搭建；
  
  

对于可以使用诸如linprog、intlinprog、quadprog等直接求解的复杂模型，现在MATLAB已经可以彻底代替Lingo实现求解了，不过在这篇文章里，还留下了很多问题没有回答，例如：规模稍大的问题，其求解效率如何？能不能得到全局最优？当约束条件变得更加复杂，例如约束之间具有因果依赖的逻辑关系时如何表示...等等。

这些问题将在这个系列的后续文章中，通过不同的代码实例一个个给予回答。

%% 1.sets
n = 2;
X = sdpvar(1,n);
%% 2.data
bulk_to_kg = [3,4];
price = [24, 16];
work_time = [12, 8];
total_milk = 50;
upper_A1 = 100;
work_time_all = 480;

%% 5.object_and_constraints
% f = X*bulk_to_kg'*price'; 向量操作，但对复杂分支结构未必通用
% c1 = sum(X)<=total_milk;

f = 0;
for i = 1 : n
    f = f + X(i)*bulk_to_kg(i)*price(i);
end
[c1, c2] = deal(0);
for i = 1 : n
    c1 = c1 + X(i);
    c2 = c2 + work_time(i)*X(i);
end
c3 = 3 * X(1);
C = [c1<=total_milk, c2<=work_time_all,c3<=upper_A1];

% solve
options = sdpsettings('solver','gurobi','verbose', 0);
diagnostics  = optimize(C,-f,options);
x = value(X)
o = -value(f)
这个代码为什么我运行出来为NaN

看看你的gurobi是不是可用的