線性規划中的單純形法與內點法（原理、步驟以及matlab實現）（四）

本篇是線性規划系列中的最后一篇，討論內點法（interior point method），相關算法在這里

原理本人也沒有搞懂，所以本文的重點在於應用

內點法不能處理等式約束，只能處理不等式約束

由對偶的相關定理我們知道如果原問題的可行解的目標函數值和對偶問題的可行解的目標函數值一致，那么該可行解是最優解

內點法的核心思想是結合原問題和對偶問題，尋找使得原問題函數值和對偶問題函數值差值等於零的可行解，通過搜尋可行域的可能解，而不僅僅是邊界

將原問題轉為karmarkar standard form：

目標：

 

步驟：

 

注：這里的兩步實現了結合原問題和對偶問題

 

注：這里的步驟作用是在結合的基礎上添加對應的約束條件

 

注：這里引入虛擬變量d，使約束條件均勻化？？？？（這里的作用不理解）

 

注：這里將所有的未知數用y表示

 

注：這里引入最后一個變量，作為目標函數（最小化），為了結合目標函數和約束條件，相應的把該變量加入每一個約束條件中，思想是系數和是0

 

步驟總結如下：

結合原問題和對偶問題——》引入相應的約束條件——》均勻化——》用一個未知數表示——》目標函數構造、對應的約束條件

下面舉例說明轉換

解答：

1.結合原問題和對偶問題

2.結合兩個問題的約束條件

 

 

3.均勻化

 

4.一個未知數表示

 

5.構造目標函數、約束條件

 

matlab實現：

由於沒有相關的內置函數所以我自己寫了相關的函數，代碼如下：

 

function y = karmarkar(f, A, b, inq, minimize, k)

%{
Input:-   1)f : The cost vector or the (row) vector containing co-eficients 
                of deceision variables in the objective function. It is a
                required parameter.
          2)A : The coefficient matrix of the left hand side of the
                constraints. it is a required parameter.
          3)b : The (row) vector containing right hand side constant of
                the constraints. It is a required parameter.
          4)inq : A (row) vector indicating the type of constraints as 1
                for >=, -1 for <= constraints. 
          5)minimize : This parameter indicates whether the objective
                function is to be minimized. minimized = 1 indicates
                a mimization problem and minimization = 0 stands for a
                maximization problem. By default it is taken as 0. It is an
                optional parameter.
%}
if minimize == 0
    f = f;
else
    f = -f;
end

for i = 1: length(inq)
    if inq(i) == 1
        A(i, :) = -A(i, :)
        b(i) = -b(i)
    end
end

%dual quetion
f_dual = b;
A_dual = A';
b_dual = f;

%change <= to = ---add s---
mn = size(A);
m = mn(1);
n = mn(2);
A = [A, eye(m)];

%change dual >= to =  ----minus s----
A_dual = [A_dual, -1 *  eye(n)];
mn = size(A);
m = mn(1);
n = mn(2);
mn_dual = size(A_dual);
m_dual = mn_dual(1);
n_dual = mn_dual(2);

A = [A, zeros(m, n)];
A_dual = [zeros(m_dual, 2 * n - n_dual), A_dual];
%add k and s
%num_pri_dual = 2 * n;
add_k_s = [ones(1, 2 * n), 1, -k];

%add constraint
num_pri = size(f);
num_dual = size(f_dual);
c_1 = [f, zeros(1, n - num_pri(2)), -1* f_dual, zeros(1, n - num_dual(2)), 0, 0];

%add d
add_d = [ones(1, 2 * n), 1, 1];

%homogenized & add object variable
A = [A, zeros(m, 1), -b', -sum([A, -b'], 2); A_dual, zeros(m_dual, 1), -b_dual', -sum([A_dual, -b_dual'], 2); 
    c_1, -sum(c_1); add_k_s, -sum(add_k_s); add_d, -sum(add_d)];
A(end, end) = 1;
disp('The Karmarkar standard form is:')
disp(A)

num_y = 2 *n + 3;
D = [1: num_y; A; ];

C = [zeros(num_y - 1, 1); 1];
E = ones(1, num_y);

 
mn =size(A);
m=mn(1);
n=mn(2);


format short e;
X1= ([E]/n)';
I=eye (n);
r=1/sqrt(n*(n-1));
alpha=(n-1)/(3*n);
X= ([E]/n)';
tol=10^(-5);
k_iteration = 0;
%while(C'*X > tol)
for i = 1:1400
	D=diag(X);
	T=A*D;
	P= [T;E];
	Q=C'*D;
	R= (I-P'/(P*P')*P)*Q';
	RN=norm (R);
	Y=X1-(r*alpha)*(R/RN);
	X=(D*Y)/(E*D*Y);
	Z=C'*X;
	k_iteration=k_iteration+1;
end
opti_x = (k+1)*X;
mn = size(f);
n = mn(2);
disp('The optimum solution is given by:')
x = opti_x(1: n, 1)
if minimize == 0
    disp(['with the maximum z = ', num2str(f*x)])
else
    disp(['with the minimum z = ', num2str(-f*x)])
end
end

 

 輸入值可看注釋

利用本函數解決本例題的代碼如下：

f = [2 1];
A = [1 1; 1 -1];
b = [5 3];
inq = [-1, -1];
karmarkar(f, A, b, inq, 0, 100)

 運行結果：

即最優解是x1=4，x2=1，z=9

下面利用linprog()檢驗：

f = [-2 -1];
A = [1 1; 1 -1];
b = [5 3];
lb = [0 0];
[x, fval] = linprog(f, A, b, [], [], lb)

 運行結果：

可以看到，最優解是一樣的