代碼 | 求解LP問題單純形法的Matlab代碼

2021-02-21 運籌學分享交流


Contents

求LP問題的基本(可行)解Matlab代碼

基本原理

Matlab代碼

初始化

是否隨機生成LP問題還是自己指定LP問題

初始化係數矩陣A和b

找到初始的基本可行解

展示線性規劃問題

求解線性規劃的單純型法

function script_LP()

求LP問題的基本(可行)解Matlab代碼基本原理

考慮如下線性規劃問題的標準型: max Z= CX s.t. {AX=b,X>=0}

假定A是m*n階矩陣,並且A的秩r(A)=m。那麼,求解該問題的最優解的的步驟: 

找到初始的基本可行解;

 判斷是否最優;

 非最優時,從一個基本可行解變換到更好的基本可行解(進基、出基、旋轉);

單純形法的本質上是每次去找更好的頂點(基本可行解),一直找到最優的頂點(基本可行解)為止。

換句話說,單純形法的要點包括以下三個部分:

第一部分,可以用之前的求解基本可行解的代碼來實現。第二和第三部分,則需要用到LP問題的典式的概念。

Matlab代碼初始化

clc
clear all
close all
format rat

是否隨機生成LP問題還是自己指定LP問題

choose = 0;

if choose

初始化係數矩陣A和b

m = 3;
n = 8;


is_opt_flag = 0;
NumMax = 10;
while is_opt_flag==0
A = randi(NumMax,m,n);
if rank(A) ==m
is_opt_flag=1;
end
end

b = randi(NumMax,m,1);
C = randi(NumMax,n,1);

[m,n] = size(A);

找到初始的基本可行解

idxs = nchoosek(1:n,m);
count = size(idxs,1);
have_no_basic_feasible = 1;
k = 1;
while k <= count && have_no_basic_feasible

BIdx = idxs(k,:);
B = A(:,BIdx);


str = [repmat('P%d,',1,m-1) 'P%d'];

if det(B) ~= 0

X = zeros(1,n);
XB = B\b;
XN = 0;
X(BIdx)=XB;


if sum(XB<0) == 0
fprintf(['\n=========================================\n1. 列向量' str '能構成基陣.\n'],BIdx);
fprintf('2. 該基陣對應的初始基本可行解:\n X= [%s].\n\n=========================================\n',rats(X));
have_no_basic_feasible = 0;
end
end
k=k+1;
end

if have_no_basic_feasible
fprintf('\n=========================================\n1. 未找到初始的基本可行解.\n');
return;
end

else
A = [1 2 1 0 0;3 2 0 1 0; 0 2 0 0 1];
b = [30;60;24];
C = [40;50;0;0;0];

[m,n] = size(A);
X = [0;0;30;60;24];
BIdx = [3,4,5];
end

展示線性規劃問題

disp('係數矩陣A是:')
disp(A)
disp('列向量b是:')
disp(b)
disp('係數向量c是:')
disp(C)

fprintf('\n=========================================\n')
disp('目標函數f=CX是:')
str = ['f=' repmat('%dx%d+',1,n-1) '%dx%d.\n'];
v =[];
for j =1:n
v =[v,C(j),j];
end
fprintf(str,v)

disp('約束組AX=b是:')
str = [repmat('%dx%d+',1,n-1) '%dx%d=%d.\n'];
for i = 1:m
v =[];
for j =1:n
v =[v, A(i,j),j];
end
fprintf(str,[v,b(i)])
end
fprintf('\n=========================================\n')

[X,Z] = simplex(A,b,C,X,BIdx);

係數矩陣A是:

1 2 1 0 0
3 2 0 1 0
0 2 0 0 1

列向量b是:
30
60
24

係數向量c是:
40
50
0
0
0


=========================================
目標函數f=CX是:
f=40x1+50x2+0x3+0x4+0x5.
約束組AX=b是:
1x1+2x2+1x3+0x4+0x5=30.
3x1+2x2+0x3+1x4+0x5=60.
0x1+2x2+0x3+0x4+1x5=24.

=========================================

function [X,Z] = simplex(A,b,C,X,BIdx)

求解線性規劃的單純型法

單純形法的基本原理,可以從不同角度去理解: 

1)從基本解的角度來理解,就是從一個基本可行解出發,去找另一個更好的基本可行解; 

2)從典式的角度來理解,就是從一個典式變換到另外一個典式; 

3)從矩陣運算的角度來理解,就是找初始基本可行解,判斷最優(檢驗數小於0),換基運算(進基、出基、旋轉運算)的過程。

[m,n] = size(A);
is_opt_flag = 0;
NIdx = setdiff(1:n,BIdx);
iter_num = 0;
while ~is_opt_flag
B = A(:,BIdx);
N = A(:,NIdx);
CB = C(BIdx);
CN = C(NIdx);

X = zeros(1,n);
XB = B\b;
XN = 0;
X(BIdx)=XB;
Z = CB'*XB;

iter_num = iter_num+1;
fprintf('第%d次單純形迭代,當前的基本可行解:\n X= [%s].\n',iter_num,rats(X));

sigma = CN'-CB'/B*N;
[max_sigma,enter_basic_idx] = max(sigma);

if max_sigma <=0
is_opt_flag = 1;
fprintf('\n=========================================\n找到該問題的最優解:\n X= [%s].\n',rats(X));
else
enter_basic_idx = NIdx(enter_basic_idx);
dk = (B\b)./(B\A(:,enter_basic_idx));

if sum(dk<=0) == m
fprintf('3. 該問題無有界解.\n');
return;
else

[d_idx,~] = find(dk>0);
[~, leave_basic_idx] = min(dk(dk>0));
leave_basic_idx = d_idx(leave_basic_idx);
leave_basic_idx = BIdx(leave_basic_idx);
BIdx = union(setdiff(BIdx,leave_basic_idx),enter_basic_idx);
NIdx = union(setdiff(NIdx,enter_basic_idx),leave_basic_idx);
end
end
end

第2次單純形迭代,當前的基本可行解:

X= [0 0  30  60  24].
第2次單純形迭代,當前的基本可行解:
X= [0  12  6  36  0].
第3次單純形迭代,當前的基本可行解:
X= [6  12  0  18  0].
第4次單純形迭代,當前的基本可行解:
X= [15 15/2  0  0  9].

=========================================
找到該問題的最優解:
X= [15  15/2  0  0  9].

Published with MATLAB® R2015b

相關焦點

  • matlab對於含nan值的數據均值計算誤解
    一直都有想把自己遇到的問題記錄下來的想法,但是一直沒有付諸行動。學習過程中,自己也慢慢發現,同樣的問題,給別人分享會讓自己對這個問題有更加深刻的認識,同時,也讓自己能夠融會貫通。從今天起,在這裡記錄自己在學習和科研過程中遇到的一些問題、生態前沿文獻分享、好用的小工具以及一些生活隨筆等等。
  • matlab數值積分1:中點規則與複合中點規則
    複合中點規則積分方法:中點規則非常粗糙,尤其是區間比較大的的時候,為了解決這個問題 matlab代碼如下:clear allf=@(x)x^2;a=0;b=2;for n = 10:10:100 h=(b-a)/(n+2); int0=0; for j = 0:n/2 xj=a+(2*j+1)*h; int0=int0+f(xj); end fprintf
  • IATA機場代碼
    如果按官方的說法,我們所說的機場三字碼叫國際航空運輸協會機場代碼,簡稱IATA代碼;四字碼叫國際民航組織機場代碼,簡稱ICAO代碼。平時工作的時候,正常人一般都叫簡稱。在IATA的三字代碼出現以前,飛行員使用的是二字代碼。大概在20世紀30年代,那時候機場還很少,但是天氣對於飛行的重要性已經被人們熟知,所以每次飛行前的必備步驟之一就是看天氣怎麼樣。
  • 幾行Matlab代碼教你上手傅立葉變換
    直接開始,如果你不熟悉Matlab,可以將代碼直接複製到編輯區,然後回車看結果就可以了。
  • 為什麼你寫的代碼連code review(代碼審查)都過不了?!
    啊,代碼啊,如流水。可是你為何又那麼美。漂亮的代碼,就像一首詩一樣,讀來讓人心動。小司機就很想做一個詩人,用代碼進行創作。寫山,寫水,寫人生。咳咳。可現實是殘酷的,小司機最近寫的代碼,反反覆覆,修修改改,可是別說優美了,連code review都過不了!!!
  • 你還敢亂寫代碼??
    代碼變壞的根源在討論什麼代碼是好代碼之前,我們先討論什麼是不好的。計算機是人造的學科,我們自己製造了很多問題,進而去思考解法。從第一行代碼開始,就是無設計的,隨意地踩著滿地的泥坑,對於旁人的眼光沒有感覺,一個人獨舞,產出的代碼,完成了需求,毀滅了接手自己代碼的人。這個就不舉例了,每個同學應該都能在自己的項目類發現這種代碼。必須形而上的思考常常,同學們聽演講,公開課,就喜歡聽一些細枝末節的'幹活'。
  • TSP問題、MTSP問題建模及求解(MATLAB)
    = zeros(1,num_brks); for kk = 1:num_brks adjust(kk) = sum(spaces == kk); end breaks = min_tour*(1:num_brks) + cumsum(adjust); end endend%代碼參考
  • 不會寫代碼的技術經理不是好CTO
    在過去的 20 年中,開發者社區呈指數級增長,開發人員的多樣性也隨之增加。換句話說,找到具有管理能力這種軟技能的開發人員其實並不困難。我非常信奉喬幫主說過的一句話:儘管技術經理不必是團隊中能力最出眾的開發人員,但至少他們應該對相關技術有所了解。當團隊成員向老闆提出技術建議時,技術經理應該能夠針對這些建議提供有價值的反饋。
  • 貨櫃尺寸、箱型及代碼對照表
    8 種常見的貨櫃及代碼1)乾貨箱:箱型代碼GP;95碼22G12)乾貨高箱:箱型代碼GH(HC/HQ);95碼25G13)掛衣箱:箱型代碼HT;95碼22V14)開頂箱:箱型代碼OT;95碼22U15)冷凍箱:箱型代碼RF;95碼22R16)冷高箱:箱型代碼RH;95碼25R17)油罐箱:箱型代碼TK;95碼22T18)框架箱:
  • 用Excel規劃求解解決購買難題
    2.分析問題,建立目標函數。3.分析問題,建立約束條件。對問題進行分析可以發現,約束條件如下:       ①商品數量為整數       ②每樣商品的數量均大於或等於1       ③總的商品價格需大於或等於2994.規劃求解:步驟1:加載規劃求解模塊:點擊文件—選項,調出"選項"對話框
  • 常用的JS 時間轉換相關代碼!
    代碼function getNumTime(num) {  let afterDate = new Date();  afterDate.setDate(afterDate.getDate() + num);  let year = afterDate.getFullYear();  let month = afterDate.getMonth
  • C語言求水仙花數代碼解析
    問題描述輸出所有的「水仙花數」,所謂的「水仙花數」是指一個三位數其各位數字的立方和等於該數本身,例如153是「水仙花數」,因為:153 = 13 + 53 + 33。問題分析根據「水仙花數」的定義,判斷一個數是否為「水仙花數」,最重要的是要把給出的三位數的個位、十位、百位分別拆分,並求其立方和(設為s),若s與給出的三位數相等, 三位數為「水仙花數」,反之,則不是。算法設計「水仙花數」是指滿足某一條件的三位數,根據這一信息可以確定整數的取值範圍是 100〜999。
  • 代碼打開workbook的兩種方法
    今天我們就講一下操作工作簿的基礎,如何用代碼打開工作簿。這一節我們介紹兩個常用的打開工作簿的方法。在VBA中打開工作薄有兩種方法:顯式打開及隱式打開。▶1.代碼及示例如下:Subworkbook_operate()    Dim wbk As Workbook    Dim fname As String' 定義工作薄對象    fname = "F:\自編小程序\轉承臺.xlsm"』注意:這裡為工作簿路徑及文件名,即fullname。
  • 如何用VBA代碼查詢兩列數據差異?
    ▼第8行至第10行代碼將表1 A列的數據存入數組aData1。第11行至第13行代碼將表2 A列的數據存入數組aData2。第14行至第17行代碼遍歷aData1的數據,作為關鍵字存入字典,並將對應的item設置為來源表的名字"表1"。第18行代碼聲明一個結果數組aRes。
  • 白話說回溯法
    一、問題引進    八皇后問題是由國際西洋棋棋手馬克斯·貝瑟爾於1848年提出的問題,是回溯算法的典型案例。  問題表述為:在8×8格的西洋棋上擺放8個皇后,使其不能互相攻擊,即任意兩個皇后都不能處於同一行、同一列或同一斜線上,問有多少種擺法。
  • 讓你徹底搞明白規劃求解!
    規劃求解加載項打勾,確定,及加載了規劃求解宏。然後根據實際問題創建基本的數據模型,比如上例中的圖片就是一個數據模型。建好數據模型後,利用規劃求解得到上例中I7和J7中數量的方法如下:1、點擊數據選項卡,最右邊找到規劃求解命令:
  • 《第一行代碼》第三版贈書來了!
    還記得在《第二行代碼》出版的時候,我也給大家送了一些書。
  • VBA常用代碼:按指定條件批量刪除Excel工作簿
    前段時間我們分享了一期小代碼,作用是一鍵提取指定文件夾下文件名,不知道大家是否還有印象:VBA常用小代碼201:批量獲取指定文件夾下文件名