matlab 线性规划 单纯形法
Posted 白水baishui
tags:
篇首语:本文由小常识网(cha138.com)小编为大家整理,主要介绍了matlab 线性规划 单纯形法相关的知识,希望对你有一定的参考价值。
先来介绍一下单纯形法,下面解释是从国科大算法最优化课程林姝老师的课件中截取的。
接下来写代码,单纯形法函数:
%% SimplexMax.m
function [x, c, z, pt, ind_B, ind_N] = SimplexMax(c, A, b, ind_B, iter_tag)
% 单纯形法求解标准形线性规划问题: max cx s.t. Ax=b x>=0
% 输入参数: c为目标函数系数, A为等式约束方程组系数矩阵, b为等式约束方程组常数项, ind_B为松弛变量索引
% 输出参数:
% x: 最优解, x中前len(ind_b)个是x_b, 后面的是x_N;
% c: 最优系数, c中从前往后分别是x1的系数、x2的系数、x3的系数......;
% z: 最优目标函数值 z_0;
% ST存储单纯形表数据;
% res_case=0表示有最优解,res_case=1表示有无界解
[m,n] = size(A); %m约束条件个数, n决策变量数
ind_N = setdiff(1:n, ind_B); %非松弛变量的索引
ST = [];
format rat
step_iter = 0;
% 循环求解
while true
x0 = zeros(n,1);
x0(ind_B) = b; %初始基可行解
cB = c(ind_B); %计算cB
Sigma = zeros(1,n);
Sigma(ind_N) = c(ind_N) - cB*A(:,ind_N); % 计算检验数
[~, k] = max(Sigma); % 选出最大检验数, 确定松弛变量索引k
Theta = b ./ A(:,k); % 计算θ
Theta(Theta<=0) = 10000;
[~, q] = min(Theta); % 选出最小θ
el = ind_B(q); % 确定出松弛变量索引el, 主元为A(q,k)
vals = [cB', ind_B', b, A, Theta];
vals = [vals; NaN, NaN, NaN, Sigma, NaN];
ST = [ST; vals];
if ~any(Sigma > 0) %此基可行解为最优解, any存在某个>0
disp('当前基可行解为最优解');
x = x0;
z = c * x;
pt = Sigma(ind_N);
res_case = 0;
return
else
disp('当前基可行解不为最优解');
x = x0;
z = c * x;
pt = Sigma(ind_N);
res_case = 0;
end % if ~any(Sigma > 0)
if all(A(:,k) <= 0) %有无界解
disp('有无界解');
x = [];
res_case = 1;
break
end % if all(A(:,k) <= 0)
step_iter = step_iter + 1;
if step_iter == iter_tag
break
end % if step_iter
% 换基
ind_B(ind_B == el) = k; % 新的松弛变量索引
ind_N = setdiff(1:n, ind_B); % 自变量索引
% 更新A和b
A(:,ind_N) = A(:,ind_B) \\ A(:,ind_N);
b = A(:,ind_B) \\ b;
A(:,ind_B) = eye(m,m);
end
end % function [x, c, z,ST,res_case] = SimplexMax(c,A,b,ind_B)
输入函数,调用单纯形法,然后输出:
%% 单纯形法求线性规划
A = [4 1 1 0; -1 1 0 1]; % 等式约束方程组(包括自变量与松弛变量)系数矩阵
b = [16; 6]; % 等式约束方程组常数项(等号右边的部分)
c = [2 3 0 0]; % 目标函数的系数项(包括自变量与松弛变量)
ind_B = [3 4]; % 松弛变量的变量索引,一般松弛变量设置地自变量更大,并紧截止自变量定义
iter_tag = 3; % 迭代次数 从1开始,n表示单纯形法运算了几轮
% SimplexMax求的是max,输出改为min的输出就行了
[x, c, z, pt, ind_B, ind_N] = SimplexMax(c, A, b, ind_B, iter_tag);
%% 求min时把c、z、pt转负数
c = -c;
z = -z;
for i=1:length(pt)
pt(i)=-pt(i);
end
%% 输出
disp('B=');
disp(A(:,ind_B));
disp('N=');
disp(A(:,ind_N));
disp('x_B=');
disp(x(ind_B));
disp('x_B的x索引=');
disp(ind_B);
disp('x_N=');
disp(x(ind_N));
disp('x_N的x索引=');
disp(ind_N);
disp('c_B的转置=');
disp(c(ind_B));
disp('c_N的转置=');
disp(c(ind_N));
disp('z_0=');
disp(z);
disp('p的转置 pT=');
disp(pt);
输出示例:
当前基可行解不为最优解
当前基可行解不为最优解
当前基可行解为最优解
B=
4 1
-1 1
N=
1 0
0 1
x_B=
2
8
x_B的x索引=
1 2
x_N=
0
0
x_N的x索引=
3 4
c_B的转置=
-2 -3
c_N的转置=
0 0
z_0=
-28
p的转置 pT=
1 2
注意出现了三行文字:
当前基可行解不为最优解
当前基可行解不为最优解
当前基可行解为最优解
这个意思是说,前两次迭代都不是最优解,可以通过控制iter_tag
来设定最大迭代次数,最后一次是最优解。
以上是关于matlab 线性规划 单纯形法的主要内容,如果未能解决你的问题,请参考以下文章
线性规划中的单纯形法与内点法(原理步骤以及matlab实现)
线性规划中的单纯形法与内点法(原理步骤以及matlab实现)