对偶单纯形法----附MATLAB函数

山·木

于 2023-12-16 10:29:14 发布

阅读量634

点赞数

文章标签： matlab 算法机器学习

本文链接：https://blog.csdn.net/m0_66511306/article/details/135010879

版权

1对偶问题的提出

1.1拉格朗日乘数法

一种解决多元函数在一定约束条件下的极值问题的方法；

将一个有 n 个变量与k个约束条件的最优化问题转换为一个有 n +k 个变量的方程组的极值问题，其变量不受任何约束；

使用方法：
在约束条件 $\varphi (x_1,x_2,x_3,...,x_n)=0$ 下，若要求 $f(x_1,x_2,x_3,...,x_n)$ 的极值，可以通过以下步骤实现。
(1)引人一个拉格朗日乘数 $\lambda$ ，作拉格朗日函数 $L(x_1,x_2,x_3,...,x_n)=f(x_1,x_2,x_3,...,x_n)+\lambda\varphi (x_1,x_2,x_3,...,x_n)$ 。

(2)此时将 $\lambda$ 视为新增变量，分别对 $x_i$ 及 $\lambda$ 求偏导，得到驻点 $x_i$ 和 $\lambda$ 。

主要思想：通过构造拉格朗日函数，将约束条件作为拉格朗日乘数的偏导数，转而求新的无约束多元函数的极值。

1.2对偶问题的生成

标准形式的线性规划问题(以下采用矩阵的形式书写)：

$z=max\mathbf{c^,x}$

$\mathbf{Ax}=\mathbf{b}$

$\mathbf{x}\geqslant 0$

引入一个新的变量 $\mathbf{y}$ ，将问题转化为：

$g=max[\mathbf{c^,x}+\mathbf{y^,}(\mathbf{b}-\mathbf{Ax})]$

满足： $\mathbf{x}\geqslant 0$

$g(\mathbf{y})=max[\mathbf{c^,x}+\mathbf{y^,}(\mathbf{b}-\mathbf{Ax})]=\mathbf{y^,b}+max[(\mathbf{c^,}-\mathbf{y^,A})\mathbf{x}]$

原问题等价于： $z=min\mathbf{y^,b}$

使得 $\mathbf{y^,A}\geqslant \mathbf{c^,}$

其转置形式为： $z=min\mathbf{b^,y}$

使得 $\mathbf{A^,y}\geqslant \mathbf{c}$

2对偶问题的性质

2.1弱对偶性

若 $x=(x_1,x_2,...,x_n)^,$ 是原问题的可行解， $y=(y_1,y_2,y_3,...,y_m)^,$ 是其对偶问题的可行解，c、b 为价值向量和资源向量，则恒有

$\mathbf{cx}\leqslant \mathbf{b^,y}$

由弱对偶性,可得出以下推论：
(1)原问题任一可行解的目标函数值是其对偶问题目标函数值的下界；反之对偶问题任一可行解的目标函数值是其原问题目标函数值的上界。
(2) 如原问题有无界解，则其对偶问题无可行解；反之对偶问题有无界解，则其原问题无可行解。
(3) 若原问题有可行解而其对偶问题无可行解，则原问题目标函数值无界；反之对偶问题有可行解而其原问题无可行解，则对偶问题的目标函数值无界。

2.2最优性

如果 $x=(x_1,x_2,...,x_n)^,$ 是原问题的可行解， $y=(y_1,y_2,y_3,...,y_m)^,$ 是其对偶问题的可行解，c、b 为价值向量和资源向量，且有

$cx=b^,y$

则 $x=(x_1,x_2,...,x_n)^,$ 是原问题的最优解， $y=(y_1,y_2,y_3,...,y_m)^,$ 是其对偶问题的最优解。

2.3强对偶性

若原问题及其对偶问题均具有可行解，则二者均具有最优解，且它们最优解的目标函数值相等。

2.4互补松驰性

在最优解中，如果对应某一约束条件的对偶变量值为非 0，则该约束条件取严格等式。反之，如果约束条件取严格不等式，则其对应的对偶变量一定为 0。

3MATLAB实现

代码如下：

function anscell = duality(A,b,c)
[RANGE10，RANGE20] = size(A);
findj0 = [];
for i = 1:1:RANGE20
    if length(find(A(:,i) == 1)) == 1&&length(find(A(:,i)))==1
        if isempty(findj0)
            ind0 = find(A(:,i) == 1);
            addfindj0 = [ind0;i];
            findj0 = [findj0 addfindj0];
        elseif~ismember(find(A(:,i) == 1),findj0(1,:))
            ind0 = find(A(:,i) == 1);
            addfindi0 = [ ind0 ;i];
            findj0 = [findj0 addfindj0];
        end
    end
end
wid = size(findj0, 2);
if wid ~= RANGE10
    error('A的格式不正确')
end
findj10 = sortrows(findj0');
base0 = findj10(:,2);
c_z = [];
cb = [];
for t = 1:1:length(base0)
    cb = [cb,c(base0(t))];
end

forj = 1:1:RANGE20
    z(j) = dot(cb',A(:,j));
    c_z=[c_z c(j)-z(j)];
end
%以上代码在simplex_method (上一篇文章)中已经使用过，这里用于检验
if max(cz)>0
    error('不满足对偶单纯形法的检验数非正条件')
end
while true
    if min(b)>=0
        printword = '最优解为';
        optans = dot(cb,b);
        ansx = zeros(1,RANGE20);
        for j = 1:1:RANGE10
            ansx(base0(j)) = b(j);
        end
        anscell = {A,b,ansx,optans,printword,base0};
        break
    end
    liindout = find(b==min(b));
    indout = liindout(1);
    if min(A(indout,:))> = 0
        printword = '原问题无可行解';
        anscell = (A, b,[],[], printword, base0)
        break
    end
    theta = inf;
    for i = 1:1:RANGE20
        if ismember(i,find(A(indout,:)< 0))
            theta = min(theta,c_z(i)/A(indout,i));
            if cz(i)/A(indout,i)== theta
                indin = i;
            end
        end
    end
    %获得 theta
    bA = [b A];
    bA(indout,:) = bA(indout,:)./bA(indout,indin + 1);
    for i = 1:1:RANGE10
        if i == indout
            continue
        end
        d = bA(i,indin + 1);
        bA(i,:) = bA(i,:)-d*(bA(indout,:));
    end

    b = bA(:,1);
    A = bA(:,2:end);
    %换入基做新基得到新的矩阵
    findj0 = [];
    for i = 1:1:RANGE20
        if length(find(A(:,i)== 1)) == 1&&length(find(A(:,i))) ==1
            if isempty(findj0)
                ind0 = find(A(:,i)== 1);
                addfindj0 = [ind0;i];
                findj0 =[findj0 addfindj0];
            elseif~ismember(find(A(:,i)==1),findj0(1,:))
                ind0 = find(A(:,i)== 1);
                addfindj0 = [ind0;i];
                findj0 = [findj0 addfindj0];
            end
        end
    end
    findjl0 = sortrows(findj0');
    base0 = findj10(:,2);
    c_z = [];
    cb = [];
    for t = 1:1:length(base0)
        cb=[cb,c(base0(t))];
    end

    for j = 1:1:RANGE20
        z(j) = dot(cb',A(:,j));
        c_z=[c_zc(j)-z(j)];
    end
    %得到新矩阵的各个量以进人下一次循环
end

参考文献《MATLAB运筹学》、《MATLAB金融风险管理》