线性规划中的单纯形法与内点法（原理、步骤以及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)

 运行结果：

可以看到，最优解是一样的