食谱问题Diet Problem

高级算法作业

1. Diet Problem的原问题及其标准型

设Corn、Milk、Bread的数量分别为x1x_1x1、x2x_2x2、x3x_3x3.
规划问题的原始形式如下：
minf=0.18x1+0.23x2+0.05x3s.t.{5000≤107x1+500x2+0x3≤500002000≤72x1+121x2+65x3≤22500≤xi≤10i=1,2,3\begin{aligned} & min \, f=0.18x_1+0.23x_2+0.05x_3 \\ &s.t. \left\{ \begin{aligned} & 5000\leq 107x_1+500x_2+0x_3 \leq 50000\\ & 2000\leq 72x_1+121x_2+65x_3 \leq 2250 \\ & 0\leq x_i \leq 10 \quad i=1,2,3\\ \end{aligned} \right. \end{aligned}minf=0.18x1+0.23x2+0.05x3s.t.⎩⎨⎧5000≤107x1+500x2+0x3≤500002000≤72x1+121x2+65x3≤22500≤xi≤10i=1,2,3
将约束拆分：
minf=0.18x1+0.23x2+0.05x3s.t.{107x1+500x2+0x3≥5000107x1+500x2+0x3≤5000072x1+121x2+65x3≥200072x1+121x2+65x3≤2250xi≤10i=1,2,3xi≥0i=1,2,3\begin{aligned} & min \, f=0.18x_1+0.23x_2+0.05x_3 \\ &s.t. \left\{ \begin{aligned} &107x_1+500x_2+0x_3 \geq 5000\\ & 107x_1+500x_2+0x_3 \leq 50000\\ & 72x_1+121x_2+65x_3 \geq 2000 \\ & 72x_1+121x_2+65x_3 \leq 2250 \\ & x_i \leq 10 \quad i=1,2,3\\ & x_i \geq 0 \quad i=1,2,3\\ \end{aligned} \right. \end{aligned}minf=0.18x1+0.23x2+0.05x3s.t.⎩⎨⎧107x1+500x2+0x3≥5000107x1+500x2+0x3≤5000072x1+121x2+65x3≥200072x1+121x2+65x3≤2250xi≤10i=1,2,3xi≥0i=1,2,3

引入松弛变量(剩余变量)x4x_4x4、x5x_5x5、x6x_6x6、x7x_7x7、x8x_8x8、x9x_9x9、x10x_{10}x10化为标准型：
minf=0.18x1+0.23x2+0.05x3+0x4+0x5+0x6+0x7+0x8+0x9+0x10s.t.{107x1+500x2+0x3−x4=5000107x1+500x2+0x3+x5=5000072x1+121x2+65x3−x6=200072x1+121x2+65x3+x7=2250x1+x8=10x2+x9=10x3+x10=10xi≥0i=1,2,...,10\begin{aligned} & min \, f=0.18x_1+0.23x_2+0.05x_3+0x_4+0x_5+0x_6+0x_7+0x_8+0x_9+0x_{10} \\ &s.t. \left\{ \begin{aligned} & 107x_1+500x_2+0x_3-x_4 =5000\\ & 107x_1+500x_2+0x_3+x_5= 50000\\ & 72x_1+121x_2+65x_3-x_6 =2000 \\ & 72x_1+121x_2+65x_3 +x_7= 2250 \\ & x_1+x_8=10 \\ & x_2+x_9=10 \\ & x_3+x_{10}=10\\ & x_i \geq 0 \quad i=1,2, ..., 10\\ \end{aligned} \right. \end{aligned}minf=0.18x1+0.23x2+0.05x3+0x4+0x5+0x6+0x7+0x8+0x9+0x10s.t.⎩⎨⎧107x1+500x2+0x3−x4=5000107x1+500x2+0x3+x5=5000072x1+121x2+65x3−x6=200072x1+121x2+65x3+x7=2250x1+x8=10x2+x9=10x3+x10=10xi≥0i=1,2,...,10

2. 对偶问题及其标准型

对上述原问题的拆分形式，其对偶问题为：
maxw=5000y1+50000y2+2000y3+2250y4+10y5+10y6+10y7s.t.{107y1+107y2+72y3+72y4+y5≤0.18500y1+500y2+121y3+121y4+y6≤0.2365y3+65y4+y7≤0.05y1,y3≥0,y2,y4,y5,y6,y7≤0\begin{aligned} & max\, w=5000y_1+50000y_2+2000y_3+2250y_4+10y_5+10y_6+10y_7 \\ &s.t. \left\{ \begin{aligned} &107y_1+107y_2+72y_3+72y_4+y_5 \leq 0.18\\ & 500y_1+500y_2+121y_3+121y_4+y_6\leq0.23 \\ & 65y_3+65y_4+y_7\leq0.05 \\ & y_1,y_3 \geq0,\, y_2,y_4, y_5,y_6,y_7\leq0\\ \end{aligned} \right. \end{aligned}maxw=5000y1+50000y2+2000y3+2250y4+10y5+10y6+10y7s.t.⎩⎨⎧107y1+107y2+72y3+72y4+y5≤0.18500y1+500y2+121y3+121y4+y6≤0.2365y3+65y4+y7≤0.05y1,y3≥0,y2,y4,y5,y6,y7≤0

引入松弛变量(剩余变量)y8y_8y8、y9y_9y9、y10y_{10}y10，对应的标准型为：
maxw=5000y1+50000y2+2000y3+2250y4+10y5+10y6+10y7+0y8+0y9+0y10s.t.{107y1+107y2+72y3+72y4+y5+y8=0.18500y1+500y2+121y3+121y4+y6+y9=0.2365y3+65y4+y7+y10=0.05y1,y3,y8,y9,y10≥0,y2,y4,y5,y6,y7≤0\begin{aligned} & max\, w=5000y_1+50000y_2+2000y_3+2250y_4+10y_5+10y_6+10y_7+0y_8+0y_9+0y_{10} \\ &s.t. \left\{ \begin{aligned} &107y_1+107y_2+72y_3+72y_4+y_5+y_8 = 0.18\\ & 500y_1+500y_2+121y_3+121y_4+y_6+y_9=0.23 \\ & 65y_3+65y_4+y_7+y_{10}=0.05 \\ & y_1,y_3,y_8,y_9,y_{10} \geq0,\, y_2,y_4, y_5,y_6,y_7\leq0\\ \end{aligned} \right. \end{aligned}maxw=5000y1+50000y2+2000y3+2250y4+10y5+10y6+10y7+0y8+0y9+0y10s.t.⎩⎨⎧107y1+107y2+72y3+72y4+y5+y8=0.18500y1+500y2+121y3+121y4+y6+y9=0.2365y3+65y4+y7+y10=0.05y1,y3,y8,y9,y10≥0,y2,y4,y5,y6,y7≤0
将变量定义域修改为yi≥0y_i\geq0yi≥0，最终，标准型为
maxw=5000y1−50000y2+2000y3−2250y4−10y5−10y6−10y7+0y8+0y9+0y10s.t.{107y1−107y2+72y3−72y4−y5+y8=0.18500y1−500y2+121y3−121y4−y6+y9=0.2365y3−65y4−y7+y10=0.05y1,y2,y3,y4,y5,y6,y7,y8,y9,y10≥0\begin{aligned} & max\, w=5000y_1-50000y_2+2000y_3-2250y_4-10y_5-10y_6-10y_7+0y_8+0y_9+0y_{10} \\ &s.t. \left\{ \begin{aligned} &107y_1-107y_2+72y_3-72y_4-y_5+y_8 = 0.18\\ & 500y_1-500y_2+121y_3-121y_4-y_6+y_9=0.23 \\ & 65y_3-65y_4-y_7+y_{10}=0.05 \\ & y_1,y_2,y_3,y_4, y_5,y_6,y_7,y_8,y_9,y_{10} \geq0\\ \end{aligned} \right. \end{aligned}maxw=5000y1−50000y2+2000y3−2250y4−10y5−10y6−10y7+0y8+0y9+0y10s.t.⎩⎨⎧107y1−107y2+72y3−72y4−y5+y8=0.18500y1−500y2+121y3−121y4−y6+y9=0.2365y3−65y4−y7+y10=0.05y1,y2,y3,y4,y5,y6,y7,y8,y9,y10≥0

3. 原问题及对偶问题的编程求解

使用MATLAB的线性规划求解工具linprog求解，代码及结果如下：

clear
clc
tic
%% 求原问题
A = [-107, -500, 0                % 系数矩阵107, 500, 0-72, -121, -6572, 121, 65];
b = [-5000, 50000, -2000, 2250];  % 常数项
lb = [0, 0, 0];      % 变量下界
ub = [10, 10, 10];   % 变量上界Aeq = [];
beq = [];
f = [0.18,0.23,0.05];
fprintf('********** 求解原问题 **********\n');
x = linprog(f,A,b,Aeq,beq,lb,ub);
fprintf('x1 = %.4f\nx2 = %.4f\nx3 = %.4f\n', x(1), x(2), x(3));
fprintf('函数值f: %.4f\n\n', x'*f');%% 求对偶问题
A = [107, 107, 72, 72, 1, 0, 0                  % 系数矩阵500, 500, 121, 121, 0, 1, 00, 0, 65, 65, 0, 0, 1];
b = [0.18, 0.23, 0.05];                         % 常数项
lb = [0, -inf, 0, -inf, -inf, -inf, -inf];      % 变量下界
ub = [inf, 0, inf, 0, 0, 0, 0];                 % 变量上界Aeq = [];
beq = [];
f = -[5000, 50000, 2000, 2250, 10, 10, 10];   % 对偶问题是最大化问题，转换为最小化问题
fprintf('********** 求解对偶问题 **********\n');
y = linprog(f,A,b,Aeq,beq,lb,ub);
fprintf('y1 = %.4f\ny2 = %.4f\ny3 = %.4f\ny4 = %.4f\ny5 = %.4f\ny6 = %.4f\ny7 = %.4f\n'..., y(1), y(2), y(3), y(4), y(5), y(6), y(7));
fprintf('函数值w: %.4f\n\n', -y'*f');
toc

可以看到，原问题与对偶问题的最优解对应的函数值相等（线性规划的强对偶性）。

4. 实验总结与心得

通过对经典线性规划问题Diet Problem的建模与求解，了解了线性规划问题的建模方式，学习了如何变换得到线性规划问题的标准型，并用求解器求解线性规划问题；
深入理解了如何将原问题转换为其对偶问题，并了解了其背后的意义（转化为另一个等价的线性规划问题）；
通过求解原问题与对偶问题，验证了线性规划的强对偶性定理，即“若原问题有最优解，那么对偶问题也有最优解，且最优值相等”。