线性规划的求解算法

线性规划的求解算法

线性规划（linear programming ）是运筹学中的一个重要分支，在现代工业、

农业、商业、交通运输、国防军事及经济管理等诸多领域都有着广泛重要的应用。在数学系的竞赛数学建模中，也多次应用线性规划来建模从而解决实际问题。在这里介绍单纯性法和对偶单纯形法两种求解线性规划的方法。

一、单纯形法算法主体思想 标准线性规划简记如下：

LP-max LP-min s.t {

Ax b x =≥ s.t {0

Ax b x =≥

这里只以LP-min 为例。

1、算法思想

单纯形法是在已知一个可行基的前提下采用的解决线性规划的算法。步骤

如下：

（1）输入初始矩阵：01020,111121,112

,1n n m m m n a a a a a a a a a +++?????

?

?

??

???

，并化为典则形式。 用R （i ）记录单位矩阵I 中元素1的位置。

（2）求{}0min

|0,1j

j a

j n t >≤≤

若t 不存在，则得到最优解；(i),1R i n x a += （i=1,2,...m ）.其他j x =0，

停。否则，转到（3）。

（3）求,1

min{|0,1}i n it it

a a i m a λ+>≤≤。 若

λ

不存在，则LP-min 无下届，所以无最优解，停；否则，求

,1min (i)|,0,1(s)i n it it a R a i m R a λ+??

=>≤≤??

??

，转到（4）。

（4）sj

sj st a a a ?，

（j=1，2....n+1）

ij

ij sj it a a a a ?-,(i=0,1,2...m;i ≠s;j=1,2,....,n+1),

(s)t R ?,转到（2）.

二、对偶单纯形法

对偶单纯形法是在已知一个正则基的条件下的求解线性规划的方法。步骤如下:

(1)输入初始矩阵：0102

0,111121,112

,1n n m m m n a a a a a a a a a +++?????

?

?

??

???

，并化为典则形式。 用R （i ）记录单位矩阵I 中元素1的位置。

(2)求{}0min

|0,1j

j a

j n s >≤≤

若s 不存在，则得到最优解；(i),1R i n x a += （i=1,2,...m ）.其他j x =0，

停。否则，转到（3）。

(3)求,1

min{|0,1}i n it it

a a i m a λ+<≤≤。 若

λ

不存在，则LP-min 无下届，所以无最优解，停；否则，求

,1min (i)|,0,1(s)i n it it a R a i m R a λ+??

=>≤≤??

??

，转到（4）。

(4)

sj

sj

st

a

a

a

?，（j=1，2....n+1）

ij ij sj it

a a a a

?-,(i=0,1,2...m;i≠s;j=1,2,....,n+1), (s)t

R?,转到（2）.

三、程序代码

1、单纯形法。

% 求解标准型线性规划:max c*x; s.t. A*x=b;x>=0

%A1是标准系数矩阵及最后一列是资源向量，C是目标函数的系数向量

% N是(初始的)基变量的下标

%M=10000 人工变量系数

% 本函数中的A是单纯形表，包括:最后一行是初始的检验数，最后一列是资源向量b %c1是基变量系数

%输出变量sol是最优解

%输出变量val是最优值，k是迭代次数

%flag1的值代表有无最优解，0无界解，1无可行解，2无穷多解，3唯一最优解

function [sol,val,k,flag1]=ssimplex(A1,C,N)

M=10000;

[mA1,nA1]=size(A1);

C1=[C,0];

val=zeros(1,length(C));

for i=1:length(N)

c1(i)=C1(N(i));

end

for i=1:nA1

a(i)=C1(i)-c1*A1(:,i);%计算初始检验数

end

A=[A1;a]; %构造初始单纯形表

[mA,nA]=size(A);

k=0; % 迭代次数

flag=1;

while flag

for i=1:(nA-1)

if A(mA,i)<=0

flag=0;

else

flag=1;

break;

end

end

if flag==0 % 已找到最优解

val1=A(1:(mA-1),nA)';

for i=1:length(N)

if (val1(i)~=0&&abs(C(N(i)))==M)

disp('无可行解');

sol=inf;val=inf;

flag3=0;

flag1=1;

break;

else

flag3=1;

end

end

if flag3

if length(find(A(mA,1:(nA-1))==0))>length(N) disp('存在无穷多最优解');

flag1=2;

else

disp('存在最优解');

flag1=3;

end

sol=c1*val1';

end

elseif flag==1

for j=1:(mA-1)

if A(j,i)<=0

flag2=1;

else

flag2=0;break;

end

end

if flag==1&&flag2==1

disp('此线性规划问题存在无界解');

sol=inf;

val=inf;

flag1=0;

flag=0; %跳出while循环

break;

end

maxq=max(A(mA,1:(nA-1)));

[m,nb]=find(A(mA,:)==maxq); %确定入基变量的纵坐标

for s=1:(mA-1)

if A(s,nb)>0

temp(s)=A(s,nA)/A(s,nb);

else

temp(s)=10000;

end

end

k=k+1;

mino=min(temp);

[n,mb]=find(temp==mino); %确定入基变量的横坐标

if length(mb)>1

mb=mb(1);

end

ab=A(mb,nb);

A2=A;

for i=1:(mA-1)

for j=1:nA

if i==mb

A(mb,j)=A2(mb,j)/ab;

else

A(i,j)=A2(i,j)-A2(i,nb)*(A2(mb,j)/ab);

end

end

end

for i=1:length(N)

if i==mb

N(i)=nb;

end

end

for i=1:length(N)

c1(i)=C(N(i));

end

for i=1:nA

A(mA,i)=C1(i)-c1*A(1:(mA-1),i);

end

end

end

if sol~=inf

for i=1:length(C)

for j=1:length(N)

if i==N(j)

val(i)=val1(j);

end

end

end

2、对偶单纯形法。

程序代码如下：

function x=lindual(c,A,b)

[n1,n2]=size(A);

A=[A,eye(n1)];c=[-c,zeros(1,n1)];

x1=[zeros(1,n2),b'];lk=[n2+1:n1+n2]; while(1)

x=x1(1:n2);

s1=[lk',b,A]

c

x1

cc=[];ci=[];

for i=1:n1

if b(i)<0

cc=[cc,b(i)];

ci=[ci,i];

end

end

nc=length(cc);

if nc==0

fprintf('达到最优解');

break

end

cliu=cc(1);

cl=ci(1);

for j=1:nc

if abs(cc(j))>abs(cliu)

cliu=cc(j);

cl=j;

end

end

cc1=[];ci1=[];

for i=1:n1+n2

if A(cl,i)<0

cc1=[cc1,A(cl,i)];

ci1=[ci1,i];

end

nc1=length(cc1);

if nc1==0

fprintf('无可行解');

break

end

cliu=c(ci1(1))/cc1(1);

cl1=ci1(1);

for j=1:nc1

if c(ci1(j))/cc1(j)

cliu=c(ci1(j))/cc1(j);

cl1=ci1(j);

end

end

b(cl)=b(cl)/A(cl,cl1);

A(cl,:)=A(cl,:)/A(cl,cl1);

for k=1:n1

if k~=cl

b(k)=b(k)-b(cl)*A(k,cl1);

A(k,:)=A(k,:)-A(cl,:).*A(k,cl1);

end

end

c=c-c(cl1).*A(cl,:);

x1(lk(cl))=0;

lk(cl)=cl1;

for kk=1:n1

x1(lk(kk))=b(kk);

end

x=x1(1:n2);

end

四、数据验证

通过课本P20页书中例题对代码进行验证如下:

s.t min212335 1232

12241 152

0,(j1,2,3,4,5) j

x x x x x x x

x x x

x x

x

+++

-+=

--+=-

-+=-

≥=

程序运行结果如下：

>> format rat

>> A=[-1 1 -1 0 0;1 2 0 -1 0;1 0 0 0 -1]; >> c=[2 1 1 0 3];

>> b=[-2;1;2];

>> x=lindual(c,A,b)

s1 =

6 -2 -1 1 -1 0 0 1 0 0

7 1 1 2 0 -1 0 0 1 0

8 2 1 0 0 0 -1 0 0 1

c =

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

x1 =

0 0 0 0 0 -2 1 2

s1 =

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

7 1 1 2 0 -1 0 0 1 0

8 2 1 0 0 0 -1 0 0 1

c =

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

x1 =

0 0 2 0 0 0 1 2

达到最优解

x = 2 0 0 1 0 经验证：结果正确！