title: 數學規划模型——線性規划問題
date: 2020-02-26 20:08:59
categories: 數學建模
tags: [MATLAB, 數學規划模型]
Matlab中線性規划的標准型
標准型
\[min \ C^{T}X \\ \ \\ s.t. \ \begin{cases} & \text{} AX<=b \quad 不等式約束\\ & \text{}Aeg*x=beg \quad 等式約束 \\ & \text{}lb<=x<=ub \quad 上下界約束(也可以當成不等式約束) \end{cases} \]
向量的內積 ,
\[c=\begin{pmatrix} C_{1}\\ C_{2}\\ ...\\ C_{n} \end{pmatrix} \quad x=\begin{pmatrix} x_{1}\\ x_{2}\\ ...\\ x_{n} \end{pmatrix},n是決策變量的個數 \]
練習題
-
min->maxm 加負號
-
不等式約束的標准是<=,>=需要轉換
-
變量如果不在約束條件，用inf與-inf巧妙轉換
Matlab 求解線性規划 的函數
[x ,fval] = linprog [ c, A, b, Aeq, beq, lb, ub, X0]
① X0 表示給定Matlab迭代求解的初始值 ( ⼀般不⽤給)
② c, A, b, Aeq, beq, lb, ub的意義和 標准型中的意義⼀致
③ 若不存在不等式約束, 可⽤ " [ ] " 替代 A和b
④ 若不存在等式約束, 可⽤ " [ ] "替代 Aeq 和 beq
⑤ 苦某個 x⽆下界或上界, 則設置lb(i)=-inf,ub(i)=+inf
⑥ 返回的 x表示⼩值處的 x取值 ; fval表示優解處時取得的最小值
7.不是所有的線性規划都有唯一解，可能無解或有無窮多的解。
8.如果求的是最大值，別忘在最后給fval加一個負號。
上⾯三個題的代碼 ：
-
[x, fval]=linprog[c, A, b, [], [], lb]
-
[x, fval]=linprog[c, A, b,Aeg, beg, lb]
-
[x, fval]=linprog[c, A, b,Aeg, beg, lb]
fval=-fval
代碼
%% Matlab求解線性規划
% [x fval] = linprog(c, A, b, Aeq, beq, lb,ub, x0)
% c是目標函數的系數向量，A是不等式約束Ax<=b的系數矩陣，b是不等式約束Ax<=b的常數項
% Aeq是等式約束Aeq x=beq的系數矩陣，beq是等式約束Aeq x=beq的常數項
% lb是X的下限，ub是X的上限，X是向量[x1,x2,...xn]' , 即決策變量。
% 迭代的初始值為x0（一般不用給）
% 更多該函數的用法說明請看講義
%% 例題1
c = [-5 -4 -6]'; % 加單引號表示轉置
% c = [-5 -4 -6]; % 寫成行向量也是可以的，不過不推薦，我們按照標准型來寫看起來比較正規
A = [1 -1 1;
3 2 4;
3 2 0];
b = [20 42 30]';
lb = [0 0 0]';
[x fval] = linprog(c, A, b, [], [], lb) % ub我們直接不寫，則意味着沒有上界的約束
% x =
% 0
% 15.0000
% 3.0000
%
% fval =
% -78
%% 例題2
c = [0.04 0.15 0.1 0.125]';
A = [-0.03 -0.3 0 -0.15;
0.14 0 0 0.07];
b = [-32 42]';
Aeq = [0.05 0 0.2 0.1];
beq = 24;
lb = [0 0 0 0]';
[x fval] = linprog(c, A, b, Aeq, beq, lb)
% x =
% 0
% 106.6667
% 120.0000
% 0
%
% fval =
% 28
% 這個題可能有多個解，即有多個x可以使得目標函數的最小值為28（不同的Matlab版本可能得到的x的值不同，但最后的最小值一定是28）
% 例如我們更改一個限定條件：令x1要大於0（注意Matlab中線性規划的標准型要求的不等式約束的符號是小於等於0）
% x1 >0 等價於 -x1 < 0，那么給定 -x1 <= -0.1 (根據實際問題可以給一個略小於0的數-0.1)，這樣能將小於號轉換為小於等於號，滿足Matlab的標准型
c = [0.04 0.15 0.1 0.125]';
A = [-0.03 -0.3 0 -0.15;
0.14 0 0 0.07
-1 0 0 0];
b = [-32 42 -0.1]';
Aeq = [0.05 0 0.2 0.1];
beq = 24;
lb = [0 0 0 0]';
[x fval] = linprog(c, A, b, Aeq, beq, lb)
% x =
% 0.1000
% 106.6567
% 119.9750
% 0
%
% fval =
% 28.0000
%% 例題3
c = [-2 -3 5]';
A = [-2 5 -1;
1 3 1];
b = [-10 12];
Aeq = ones(1,3);
beq = 7;
lb = zeros(3,1);
[x fval] = linprog(c, A, b, Aeq, beq, lb)
fval = -fval % 注意這個fval要取負號（原來是求最大值，我們添加負號變成了最小值問題）
% x =
% 6.4286
% 0.5714
% 0
% fval =
% -14.5714
% fval =
% 14.5714
%% 多個解的情況
% 例如 ： min z = x1 + x2 s.t. x1 + x2 >= 10
c = [1 1]';
A = [-1 -1];
b = -10;
[x fval] = linprog(c, A, b) % Aeq, beq, lb和ub我們都沒寫，意味着沒有等式約束和上下界約束
% x有多個解時，Matlab會給我們返回其中的一個解
%% 不存在解的情況
% 例如 ： min z = x1 + x2 s.t. x1 + x2 = 10 、 x1 + 2*x2 <= 8、 x1 >=0 ，x2 >=0
c = [1 1]';
A = [1 2];
b = 8;
Aeq = [1 1];
beq = 10;
lb = [0 0]';
[x fval] = linprog(c, A, b, Aeq, beq, lb) % Linprog stopped because no point satisfies the constraints.（沒有任何一個點滿足約束條件）
線性規划的典型例題
例題1
設備有效台：
你有100台設備，每月工作22天，每天工作8小時，設備每天只有5小時有效運作（生產），則你的設備有效台時=5×22×100=11000台時，你的工作台時=8×22×100=17600台時，設備有效作業率=5/8=62.5%
%% 生產決策問題
format long g %可以將Matlab的計算結果顯示為一般的長數字格式（默認會保留四位小數，或使用科學計數法）
% (1) 系數向量
c = zeros(9,1); % 初始化目標函數的系數向量全為0
c(1) = 1.25 -0.25 -300/6000*5; % x1前面的系數是c1
c(2) = 1.25 -0.25 -321/10000*7;
c(3) = -250 / 4000 * 6;
c(4) = -783/7000*4;
c(5) = -200/4000 * 7;
c(6) = -300/6000*10;
c(7) = -321 / 10000 * 9;
c(8) = 2-0.35-250/4000*8;
c(9) = 2.8-0.5-321/10000*12-783/7000*11;
c = -c; % 我們求的是最大值，所以這里需要改變符號
% (2) 不等式約束
A = zeros(5,9);
A(1,1) = 5; A(1,6) = 10;
A(2,2) = 7; A(2,7) = 9; A(2,9) = 12;
A(3,3) = 6; A(3,8) = 8;
A(4,4) = 4; A(4,9) = 11;
A(5,5) = 7;
b = [6000 10000 4000 7000 4000]';
% (3) 等式約束
Aeq = [1 1 -1 -1 -1 0 0 0 0;
0 0 0 0 0 1 1 -1 0];
beq = [0 0]';
%（4）上下界
lb = zeros(9,1);
% 進行求解
[x fval] = linprog(c, A, b, Aeq, beq, lb)
fval = -fval
% fval =
% 1146.56650246305
% 注意，本題應該是一個整數規划的例子，我們在后面的整數規划部分再來重新求解。
intcon = 1:9;
[x,fval]=intlinprog(c,intcon,A,b,Aeq,beq,lb)
fval = -fval
例題2
%% 投料問題
clear,clc
format long g %可以將Matlab的計算結果顯示為一般的長數字格式（默認會保留四位小數，或使用科學計數法）
% (1) 系數向量
a=[1.25 8.75 0.5 5.75 3 7.25]; % 工地的橫坐標
b=[1.25 0.75 4.75 5 6.5 7.25]; % 工地的縱坐標
x = [5 2]; % 料場的橫坐標
y = [1 7]; % 料場的縱坐標
c = []; % 初始化用來保存工地和料場距離的向量 (這個向量就是我們的系數向量）
for j =1:2
for i = 1:6
c = [c; sqrt( (a(i)-x(j))^2 + (b(i)-y(j))^2)]; % 每循環一次就在c的末尾插入新的元素
end
end
% (2) 不等式約束
A =zeros(2,12);
A(1,1:6) = 1;
A(2,7:12) = 1;
b = [20,20]';
% (3) 等式約束
Aeq = zeros(6,12);
for i = 1:6
Aeq(i,i) = 1; Aeq(i,i+6) = 1;
end
% Aeq = [eye(6),eye(6)] % 兩個單位矩陣橫着拼起來
beq = [3 5 4 7 6 11]'; % 每個工地的日需求量
%（4）上下界
lb = zeros(12,1);
% 進行求解
[x fval] = linprog(c, A, b, Aeq, beq, lb)
x = reshape(x,6,2) % 將x變為6行2列便於觀察（reshape函數是按照列的順序進行轉換的，也就是第一列讀完，讀第二列，即x1對應x_1,1，x2對應x_2,1）
% fval =
% 135.281541790676