数学规划模型 - Starslayerx/MathModeling GitHub Wiki
由于规划问题使用lingo,这里为相关文档
-
算术运算符
^ 幂 x 乘 / 除 + 加 - 减 优先级为:-(取反) ^ x / + -
-
逻辑运算符
#not# !取反 #eq# == #ne# != #gt# > #ge# >= #lt# < #le# <= #and# and #or# or 优先级为:not eq ne... le and or
字母缩写:g:greater, l:less, e:equal, t:than, n:not -
关系运算符
<= = >=
Lingo并不支持严格小于:A < B,不过可以使用A+e <= B,其中e为一个很小的正数.
-
数学函数
数学函数 解释 @abs(x) 绝对值 @sin(x) sin @cos(x) cos @tan(x) tan @exp(x) 返回e的x次方 @log(x) 返回x的自然对数 @lgm(x) 返回x的gamma函数的自然对数 @sign(x) 若x<0,返回-1,否则,返回1 @floor(x) 返回x整数部分 @smax(x1,x2,...,xn) 返回最大值 @smin(x1,x2,...,xn) 返回最小值 -
金融函数
@FPA(I, N)
单位时间利率为I,从下个时刻开始连续N个时间段支付。每个时段单位费用 -
概率函数
概率函数 解释 @pbn(p,n,x) 二项分布的累积分布函数,n或x不为整数时,使用线性插值法计算 @pcx(n,x) 自由度为n的X^2分布的累积分布函数 @peb(a,x) 当符合为a,服务器系统有x个服务器且允许无穷排队时的Erlang繁忙概率 @pel(a,x) 当符合为a,服务器系统有x个服务器且不允许排队时的Erlang繁忙概率 @pfd(n,d,x) 自由度为n和d的F分布的累积分布函数 @pfs(a,x,c) 当符合上限为a,顾客数为c,平行服务器数量为x时,有限源的Poisson服务器系统的等待或返修顾客数的期望值。a为顾客数乘以平均服务时间,再除以返修时间。当c或x不为整数时使用线性插值法机算 @phg(pop,g,n,x) 超几何分布的累积分布函数.pop为产品总数,g为正品数,从中任意取出n件。都可以为非整数,采用线性插值计算 @ppl(a,x) Poisson分布的线性损失函数,即返回max(0,z-x)的期望值,随机变量x服从均值为a的Poisson分布 @pps(a,x) 均值为a的Poisson分布的累积分布函数,当x不为整数时使用线性插值法计算 @psl(x) 单位正态线性损失函数,返回max(0,z-x)的期望值,x服从标准正态分布 @psn(x) 标准正态分布的累积分布函数 @ptd(n,x) 自由度为n的t分布的累积分布函数 @qrand(seed) 产生服从(0,1)区间的随机数 @rand(seed) 返回0~1之间的伪随机数: U(I+1) =@rand[U(I)] -
变量界定函数
对变量范围进行界定界定函数 解释 @bin(x) 0或1 @bnd(L,x,U) L<=x<=U @free(x) x可以为任意实数(Lingo默认x>=0) @gin(x) 限制x为整数 -
集合操作函数
函数 功能 @in() 元素是否在指定集合中 @index([a], b) 返回a在b集合中的索引,找不到会报错 @wrap(index, limit) 返回index (mod) limit @size(set_name) 返回集合成员个数 -
集合循环函数
@function(set_operator(set_name |condition: expression))- set_operator: 集合函数名
- set_name: 数据集合名
- |condition: 条件
- expressino: 表达式
常见集合函数 解释 @for(set_name: constraint_expressions) 对集合中每个元素独立地生成约束,由约束表达式概述 @max(set_name: expression) 返回集合上表达式的最大值 @min(set_name: expression) 返回集合上表达式的最小值 @sum(set_name: expression) 返回集合表达式和和 @size(set_name) 返回元素个数 @in(set_name, set_element) 若set_name中包含set_element则返回1,否则返回0 -
输入输出函数
函数 解释 @file('filename') 读取文件 @text() 输出到文件 @ole() 从Excel中引入或输出数据的接口函数 @status() 返回Lingo求解模型后状态 @dual() 返回变量的判别数或约束行的影子价格 Lingo 求解的状态有
0 GlobalOptimum (全局最优解)
1 Infeasible (不可行)
2 Unbounded (无界)
3 Undetermined (不确定)
4 Feasible (可行)
5 InfeasibleorUnbounded (通常需要关闭“预处理”选项后重新求解模型)
6 LocalOptimum (局部最优解)
7 LocallyInfeasible (局部不可行,可能有解,但Lingo没有找到)
8 Cutoff (目标函数的截断值被达到)
9 NumericError (求解器在某约束中遇到无定义的算术运算而停止)
- 产品信息
每桶牛奶甲(A1) 乙(A2) 时间 12h 产量 3kg 价格 24¥/kg
可得: 利润 = 324x1 + 416x2
即 z = 72x1 + 64x2
-
约束
变量 值 时间 最大480h 甲产量 10kg 牛奶 最大50桶
即
x1 + x2 <= 50 (milk)
12x1 + 8x2 <= 480 (time)
3x1 <= 100 (甲产能)
x1 >=0, x2 >= 0
model:
[obj] max = 72*x1+64*x2;
[milk] x1+x2<50;
[time] 12*x1+8*x2<480;
[cpct] 3*x1<100;
end
- model: 表示模型开始
- end: 模型结束
- []: 括号内的为注释,对这一行的说明
运行结果如下:
Global optimal solution found.(全局最优解)
Objective value: 3360.000 (目标函数值)
Infeasibilities: 0.000000
Total solver iterations: 2
Elapsed runtime seconds: 0.07
Model Class: LP
Total variables: 2
Nonlinear variables: 0
Integer variables: 0
Total constraints: 4
Nonlinear constraints: 0
Total nonzeros: 7
Nonlinear nonzeros: 0
Variable Value Reduced Cost
X1 20.00000 0.000000
X2 30.00000 0.000000
Row Slack or Surplus Dual Price
(资源剩余量) (最优解下资源增加1目标函数的变化,影子价格)
OBJ 3360.000 1.000000
MILK 0.000000 48.00000
TIME 0.000000 2.000000
CPCT 40.00000 0.000000
- 资源剩余0为 紧约束
- 增加紧约束的资源会增加利润 影子价格
敏感度分析运行结果:
Ranges in which the basis is unchanged:
Objective Coefficient Ranges:
Current Allowable Allowable
Variable Coefficient Increase Decrease
(当前系数) (最优解不变的情况下,系数的变化范围)
X1 72.00000 24.00000 8.000000
X2 64.00000 8.000000 16.00000
Righthand Side Ranges:
Current Allowable Allowable
Row RHS Increase Decrease
(当前右端项,总资源) (影子价格有意义的范围,只是范围内一定有意义,超出范围后并不确定)
MILK 50.00000 10.00000 6.666667
TIME 480.0000 53.33333 80.00000
CPCT 100.0000 INFINITY 40.00000
例1解答:
- 每增加1桶牛奶,利润增加48¥,所以应该投资,最多买10桶
- 最多每小时不超过2¥,但最多增加53.3h
- 30*3 = 90; x1(max) = 72+24 =96 > 90 在范围内,不用改变生产计划
产品深加工
| A1->B1 | A2->B2 | |
|---|---|---|
| 时间 | 2h | 3h |
| 产量 | 0.8kg | 0.75kg |
| 利润 | 44¥/kg | 32¥/kg |
假设: x1 kg A1, x2 kg A2, x3 kgB1, x4 kg B2, x5 kg A1->B1, x6 kg A2->B2 (x1~x2为最后得到的A产品产量, x5~x6用来简化问题)
[注意]: 这次设x为kg,上次x为桶数
目标函数: z = 24x1 + 16x2 + 44x3 + 32x4 - 3x5 - 3x6
-
约束
原料供应 (x1+x5) / 3 + (x2+x6) / 4 <= 50 劳动时间 4(x1+x5) + 2(x2+x6) + 2x5 + 2x6 <= 480 设备产能 x1 + x5 <= 100 附加约束 x3 = 0.8x5, x4 = 0.75x6 非负约束 x1~x6 >= 0
lingo程序
[注意]: 约束条件式子单位意义要搞清楚,例如[milk]式子由于两边乘了12,单位为(1/12 桶),灵敏度分析时的milk单位不要弄错了
model:
[obj] max = 24*x1 + 16*x2 + 44*x3 + 32*x4 - 3*x5 - 3*x6;
[milk] (x1+x5)*4 + (x2 + x6)*3 <= 50*12;
[time] 4*(x1+x5) + 2*(x2+x6) + 2*x5 + 2*x6 <= 480;
[cpct] x1 + x5 <= 100;
x3 = 0.8*x5;
x4 = 0.75*x6;
end
求解(含敏感度分析)
Global optimal solution found.
Objective value: 3460.800
Infeasibilities: 0.000000
Total solver iterations: 2
Elapsed runtime seconds: 0.03
Model Class: LP
Total variables: 6
Nonlinear variables: 0
Integer variables: 0
Total constraints: 6
Nonlinear constraints: 0
Total nonzeros: 20
Nonlinear nonzeros: 0
Variable Value Reduced Cost
X1 0.000000 1.680000
X2 168.0000 0.000000
X3 19.20000 0.000000
X4 0.000000 0.000000
X5 24.00000 0.000000
X6 0.000000 1.520000
Row Slack or Surplus Dual Price
OBJ 3460.800 1.000000
MILK 0.000000 3.160000
TIME 0.000000 3.260000
CPCT 76.00000 0.000000
5 0.000000 44.00000
6 0.000000 32.00000
Ranges in which the basis is unchanged:
Objective Coefficient Ranges:
Current Allowable Allowable
Variable Coefficient Increase Decrease
X1 24.00000 1.680000 INFINITY
X2 16.00000 8.150000 2.100000
X3 44.00000 19.75000 3.166667
X4 32.00000 2.026667 INFINITY
X5 -3.000000 15.80000 2.533333
X6 -3.000000 1.520000 INFINITY
Righthand Side Ranges:
Current Allowable Allowable
Row RHS Increase Decrease
MILK 600.0000 120.0000 280.0000
TIME 480.0000 253.3333 80.00000
CPCT 100.0000 INFINITY 76.00000
5 0.000000 INFINITY 19.20000
6 0.000000 INFINITY 0.000000
生产计划见结果,不再阐述
- 增加每桶牛奶利润: 12*3.16 = 37.92¥, 增加每小时工作时间利润: 3.26¥,净利润分别为: 7.92¥和0.26¥,因此该尽量多的增加牛奶数量,牛奶最多可以增加120 * 1/12 = 10桶,所以我们用150元买5桶30元的牛奶利润最高,利润为5*37.92 = 189.6¥,净利润为5*7.92 = 39.6¥。
- 根据x3 x4的敏感度分析知,B1下降10% 或 B2上升10%生产计划不一定为最优的。此时,若要继续分析,需要将x3(B1)的系数改成44 - 4.4 = 39.6¥,或将x4(B2)系数改成32 - 4.4 = 27.6¥,重新计算该模型,根据模型调整生产计划
- 根据x1的Reduce Cost严格大于0,首先:最优解中x1取值一定为0,其次:如果限定x1取值大于等于某个正数,则x1从0开始增加一个单位时,目标函数必将减少1.68。因此,10kg的A1合同,必须满足,会使公司利润减少1.68 * 10 = 16.8元。(Reduce Cost: 最优单纯形表中判别数所在行的变量的系数,表示当变量有微小变动时,目标函数的变化率。)
例2的程序:
!定义集合;
sets:
cang/1..3/:WET,VOL;
wu/1..4/:w,v,p;
link(wu,cang):x;
endsets
!给集合赋值;
data:
WET = 10, 16, 8;
VOL = 6800, 8700, 5300;
w = 18, 15, 23, 12;
v = 480, 650, 580, 390;
p = 3100, 3800, 3500, 2850;
enddata
max = @sum(wu(i):p(i) * @sum(cang(j):x(i,j)));
@for(wu(i):@sum(cang(j):x(i,j))<w(i));
@for(cang(j):@sum(wu(i):x(i,j))<WET(j));
@for(cang(j):@sum(wu(i):v(i)*x(i,j))<VOL(j));
@for(cang(j):
@for(cang(k)|k #gt# j:
@sum(wu(i):x(i,j)/WET(j))=@sum(wu(i):x(i,k)/WET(k)));
);
end
输出如下:
Global optimal solution found.
Objective value: 121515.8
Infeasibilities: 0.000000
Total solver iterations: 12
Elapsed runtime seconds: 0.14
Model Class: LP
Total variables: 12
Nonlinear variables: 0
Integer variables: 0
Total constraints: 14
Nonlinear constraints: 0
Total nonzeros: 72
Nonlinear nonzeros: 0
Variable Value Reduced Cost
WET( 1) 10.00000 0.000000
WET( 2) 16.00000 0.000000
WET( 3) 8.000000 0.000000
VOL( 1) 6800.000 0.000000
VOL( 2) 8700.000 0.000000
VOL( 3) 5300.000 0.000000
W( 1) 18.00000 0.000000
W( 2) 15.00000 0.000000
W( 3) 23.00000 0.000000
W( 4) 12.00000 0.000000
V( 1) 480.0000 0.000000
V( 2) 650.0000 0.000000
V( 3) 580.0000 0.000000
V( 4) 390.0000 0.000000
P( 1) 3100.000 0.000000
P( 2) 3800.000 0.000000
P( 3) 3500.000 0.000000
P( 4) 2850.000 0.000000
X( 1, 1) 0.000000 400.0000
X( 1, 2) 0.000000 57.89474
X( 1, 3) 0.000000 400.0000
X( 2, 1) 7.000000 0.000000
X( 2, 2) 0.000000 239.4737
X( 2, 3) 8.000000 0.000000
X( 3, 1) 3.000000 0.000000
X( 3, 2) 12.94737 0.000000
X( 3, 3) 0.000000 0.000000
X( 4, 1) 0.000000 650.0000
X( 4, 2) 3.052632 0.000000
X( 4, 3) 0.000000 650.0000
Row Slack or Surplus Dual Price
1 121515.8 1.000000
2 18.00000 0.000000
3 0.000000 300.0000
4 7.052632 0.000000
5 8.947368 0.000000
6 0.000000 3500.000
7 0.000000 3265.789
8 0.000000 0.000000
9 510.0000 0.000000
10 0.000000 3.421053
11 100.0000 0.000000
12 0.000000 0.000000
13 0.000000 0.000000
14 0.000000 -28000.00
类似的指派问题(竞赛指南Page59):

model:
sets:
supply/1,2/:s;
demand/1,2,3/:d;
link(supply, demand):road, sd;
endsets
data:
road = 10, 5, 6, 4, 8, 12;
d = 50, 70, 40;
s = 60, 100;
enddata
[obj]min = @sum(link(i,j): road(i,j)*sd(i,j));
@for(demand(j):@sum(supply(i):sd(i,j)) = d(j));
@for(supply(i):@sum(demand(j):sd(i,j)) < s(i));
end
该问题为整数规划问题(Intger Programming, IP),在lingo中使用@gin(x)将变量约束为整数
model:
[obj] max = 2*x1 + 3*x2 + 4*x3;
[material] 1.5*x1 + 3*x2 + 5*x3 <= 600;
[time] 280*x1 + 250*x2 + 400*x3 <= 60000;
@gin(x1);
@gin(x2);
@gin(x3);
end
得到生产车辆类型的最优解
对于第二问可以使用非线性规划(Non-Linear Programming)的方法求解
x1 * (x1 - 80) >= 0
x2 * (x2 - 80) >= 0
x3 * (x3 - 80) >= 0
将以上约束添加入程序中:
model:
[obj] max = 2*x1 + 3*x2 + 4*x3;
[material] 1.5*x1 + 3*x2 + 5*x3 <= 600;
[time] 280*x1 + 250*x2 + 400*x3 <= 60000;
x1*(x1-80) >= 0;
x2*(x2-80) >= 0;
x3*(x3-80) >= 0;
@gin(x1);
@gin(x2);
@gin(x3);
end
一般非线性规划求解非常困难,还可以将问题所用情况列举出来然后分解求解模型的方法来得到最优解;此外,还可以引入0,1变量建立整数规划模型
80*y1 <= x1 <= M*y1, @bin(y1);
80*y2 <= x2 <= M*y2, @bin(y2);
80*y3 <= x3 <= M*y3, @bin(y3);
(M为相当大的整数)
例2 原油采购加工
该题的例如计算涉及到一个分段函数
- 可以将公式改为
c(x) = 10*x1 + 8*x2 + 6*x3, x = x1 + x2 + x3
相应的有如下约束,当x1 < 500, x2 =0, 当x1 = 500, 0 < x2 <= 500 (x2与x3同理),由此得到如下公式
(x1 - 500)*x2 = 0
(x2 - 500)*x3 = 0 0 <= x1,x2,x3 <= 500
lingo程序如下
model:
[obj] max = 4.8*x11 + 4.8*x21 + 5.6*x12 + 5.6*x22 - 10*x1 - 8*x2 - 6*x3;
x11 + x12 < x + 500;
x21 + x22 < 1000;
0.5 * x11 - 0.5 * x21 > 0;
0.4 * x12 - 0.6 * x22 > 0;
x = x1 + x2 + x3;
(x1 - 500) * x2 = 0;
(x2 - 500) * x3 = 0;
x1 < 500;
x2 < 500;
x3 < 500;
end
但这次求得的是局部最优解
Local optimal solution found.
Objective value: 4800.000
Infeasibilities: 0.000000
Total solver iterations: 8
Elapsed runtime seconds: 0.14
Model Class: QP
Total variables: 8
Nonlinear variables: 3
Integer variables: 0
Total constraints: 11
Nonlinear constraints: 2
Total nonzeros: 27
Nonlinear nonzeros: 2
Variable Value Reduced Cost
X11 500.0000 0.000000
X21 500.0000 0.000000
X12 0.000000 0.000000
X22 0.000000 0.4000000
X1 0.000000 0.4000000
X2 0.000000 0.000000
X3 0.000000 0.000000
X 0.000000 0.000000
Row Slack or Surplus Dual Price
OBJ 4800.000 1.000000
2 0.000000 9.600000
3 500.0000 0.000000
4 0.000000 -9.600000
5 0.000000 -10.00000
6 0.000000 9.600000
7 0.000000 -0.3200000E-02
8 0.000000 -0.7200000E-02
9 500.0000 0.000000
10 500.0000 0.000000
11 500.0000 0.000000
要自行修改使用全局solver,之后的求解如下
Global optimal solution found.
Objective value: 5000.002
Objective bound: 5000.002
Infeasibilities: 0.7253991E-06
Extended solver steps: 3
Total solver iterations: 116
Elapsed runtime seconds: 0.19
Model Class: QP
Total variables: 8
Nonlinear variables: 3
Integer variables: 0
Total constraints: 11
Nonlinear constraints: 2
Total nonzeros: 27
Nonlinear nonzeros: 2
Variable Value Reduced Cost
X11 0.000000 0.9000000
X21 0.000000 0.000000
X12 1500.000 0.000000
X22 1000.000 0.000000
X1 500.0000 0.000000
X2 499.9991 0.000000
X3 0.8517036E-03 0.000000
X 1000.000 0.000000
Row Slack or Surplus Dual Price
OBJ 5000.002 1.000000
2 0.000000 7.000000
3 0.000000 3.500000
4 0.000000 -2.600000
5 0.000000 -3.500000
6 0.000000 7.000000
7 0.000000 -0.6000010E-02
8 -0.7253991E-06 -1174.117
9 0.000000 0.000000
10 0.8517036E-03 0.000000
11 499.9991 0.000000
可见局部最优解并不是全局最优解,此外也可以将模型转化为线性约束(该方法更容易求解)

500 * y2 <= x1 <= 500 * y1
500 * y3 <= x2 <= 500 * y2
x1 <= 500 * y3
y1, y2, y3 = 0 或 1
还可以直接处理该分段函数,下面是python绘制程序
x1 = np.linspace(1, 500, 1000)
y1 = np.linspace(1, 5000, 1000)
x2 = np.linspace(500, 1000, 1000)
y2 = np.linspace(5000, 9000, 1000)
x3 = np.linspace(1000, 1500, 1000)
y3 = np.linspace(9000, 12000, 1000)
fig = plt.figure()
# 绘图主题
# plt.style.use('seaborn-whitegrid')
# 添加刻度
plt.xticks([0, 500, 1000, 1500])
plt.yticks([5000, 9000, 12000], [5000, 9000, 12000])
# 控制坐标范围(这样坐标轴不会有多出来的部分,可以去掉比较一下)
plt.xlim(0, 1500)
plt.ylim(0, 12000)
# 添加坐标轴
plt.xlabel('$x$', fontproperties='SimHei')
plt.ylabel('$c(x)$', fontproperties='SimHei')
# 绘制分段函数
plt.plot(x1, y1)
plt.plot(x2, y2)
plt.plot(x3, y3)
# 绘制虚线
plt.vlines(500, 0, 5000, colors = "grey", linestyles = "dashed")
plt.vlines(1000, 0, 9000, colors = "grey", linestyles = "dashed")
# plt.vlines(1500, 0, 12000, colors = "grey", linestyles = "dashed")
plt.hlines(5000, 0, 500, colors= 'grey', linestyles = 'dashed')
plt.hlines(9000, 0, 1000, colors= 'grey', linestyles = 'dashed')
# plt.hlines(12000, 0, 1500, colors= 'grey', linestyles = 'dashed')
plt.show()
下面为该问题的解决方法:

该类问题关键就是处理分段函数,最后一种方法更具一般性,通用步骤如下

例1 混合泳接力队的选拔(指派问题 Assignment)
这里使用0-1规划,将X为0或1,表示该人是否参加,c为队员i的j类游泳姿势成绩
model:
sets:
person/1..5/;
position/1..4/;
link(person, position): c,x;
endsets
data:
c = 66.8, 57.2, 78, 70, 67.4,
75.6, 66, 67.8, 74.2, 71,
87, 66.4, 84.6, 69.6, 83.8,
58.6, 53, 59.4, 57.2, 62.4;
enddata
min = @sum(link:c*x);
@for(person(i):@sum(position(j):x(i,j)) <= 1;);
@for(position(i):@sum(person(j):x(j,i)) = 1;);
@for(link:@bin(x));
end
根据输出即得解
例2 选课策略
- 0-1规划中,一般使用0代表不选中,而1代表选中(选择策略),其次0和1也可以组成逻辑关系,例如,要有B必须先有A,那么可以使用
B<=A来表示其逻辑关系 - 多目标规划中,可以使用加权的方法得到一个新的单目标规划
- 选择关系可以使用
w1 + w2 + ... + wn = 1和wi = {0,1}来表示,即wi中一个为1其余都为0
**例3 **