MATLAB和Python解线性规划

Canethui 2020-03-03

MATLAB和Python解线性规划

  1. 作可行域的方法
//画可行域的方法
[X,Y]=meshgrid(0:0.1:100,0:0.1:100);
idx=(X+Y>=10)&(-2*X+2*Y<=10)&(-4*X+2*Y<=20)&(X+4*Y>=20);
x=X(idx);
y=Y(idx);
k=convhull(x,y);
fill(x(k),y(k),'c');

//m文件,自定义函数的保存路径
G:\MATLAB\toolbox\shared\maputils
画完图后再调用xyplot函数
  1. 使y轴放置在中间, 也就是上面的xyplot函数

    %作用:将Y坐标轴放在中间
    function xyplot(x,y)
    % PLOT
    if nargin>0
    if nargin == 2
    plot(x,y);
    else
    display(' Not 2D Data set !')
    end
    end
    hold on;
    % GET TICKS
    X=get(gca,'Xtick');
    Y=get(gca,'Ytick');
    % GET LABELS
    XL=get(gca,'XtickLabel');
    YL=get(gca,'YtickLabel');
    % GET OFFSETS
    Xoff=diff(get(gca,'XLim'))./40;
    Yoff=diff(get(gca,'YLim'))./40;
    % DRAW AXIS LINEs
    plot(get(gca,'XLim'),[0 0],'k');
    plot([0 0],get(gca,'YLim'),'k');
    % Plot new ticks
    for i=1:length(X)
    plot([X(i) X(i)],[0 Yoff],'-k');
    end;
    for i=1:length(Y)
    plot([Xoff, 0],[Y(i) Y(i)],'-k');
    end;
    % ADD LABELS
    text(X,zeros(size(X))-2.*Yoff,XL);
    text(zeros(size(Y))-3.*Xoff,Y,YL);
    box off;
    % axis square;
    axis off;
    set(gcf,'color','w');
    set(gca,'FontSize',20);
  2. MATLAB直接求解线性规划

    >> f=[2,-1];
    >> A=[-1,-1;-2,2;-4,2;-1,-4];
    >> b=[-10;10;20;-20];
    >> lb=zeros(2,1);
    >> [x,fval]=linprog(f,A,b,[],[],lb,[])
  3. Python代码实现单纯形法

    # coding=utf-8
    # 单纯形法的实现,只支持最简单的实现方法
    # 且我们假设约束矩阵A的最后m列是可逆的
    # 这样就必须满足A是行满秩的(m*n的矩阵)
    
    import numpy as np
    
    
    class Simplex(object):
        def __init__(self, c, A, b):
            # 形式 minf(x)=c.Tx
            # s.t. Ax=b
            self.c = c
            self.A = A
            self.b = b
    
        def run(self):
            c_shape = self.c.shape
            A_shape = self.A.shape
            b_shape = self.b.shape
            assert c_shape[0] == A_shape[1], "Not Aligned A with C shape"
            assert b_shape[0] == A_shape[0], "Not Aligned A with b shape"
    
            # 找到初始的B,N等值
            end_index = A_shape[1] - A_shape[0]
            N = self.A[:, 0:end_index]
            N_columns = np.arange(0, end_index)
            c_N = self.c[N_columns, :]
            # 第一个B必须是可逆的矩阵,其实这里应该用算法寻找,但此处省略
            B = self.A[:, end_index:]
            B_columns = np.arange(end_index, A_shape[1])
            c_B = self.c[B_columns, :]
    
            steps = 0
            while True:
                steps += 1
                print("Steps is {}".format(steps))
                is_optim, B_columns, N_columns = self.main_simplex(B, N, c_B, c_N, self.b, B_columns, N_columns)
                if is_optim:
                    break
                else:
                    B = self.A[:, B_columns]
                    N = self.A[:, N_columns]
                    c_B = self.c[B_columns, :]
                    c_N = self.c[N_columns, :]
    
        def main_simplex(self, B, N, c_B, c_N, b, B_columns, N_columns):
            B_inverse = np.linalg.inv(B)
            P = (c_N.T - np.matmul(np.matmul(c_B.T, B_inverse), N)).flatten()
            if P.min() >= 0:
                is_optim = True
                print("Reach Optimization.")
                print("B_columns is {}".format(B_columns))
                print("N_columns is {}".format(sorted(N_columns)))
                best_solution_point = np.matmul(B_inverse, b)
                print("Best Solution Point is {}".format(best_solution_point.flatten()))
                print("Best Value is {}".format(np.matmul(c_B.T, best_solution_point).flatten()[0]))
                print("\n")
                return is_optim, B_columns, N_columns
            else:
                # 入基
                N_i_in = np.argmin(P)
                N_i = N[:, N_i_in].reshape(-1, 1)
                # By=Ni, 求出基
                y = np.matmul(B_inverse, N_i)
                x_B = np.matmul(B_inverse, b)
                N_i_out = self.find_out_base(y, x_B)
                tmp = N_columns[N_i_in]
                N_columns[N_i_in] = B_columns[N_i_out]
                B_columns[N_i_out] = tmp
                is_optim = False
    
                print("Not Reach Optimization")
                print("In Base is {}".format(tmp))
                print("Out Base is {}".format(N_columns[N_i_in]))  # 此时已经被换过去了
                print("B_columns is {}".format(sorted(B_columns)))
                print("N_columns is {}".format(sorted(N_columns)))
                print("\n")
                return is_optim, B_columns, N_columns
    
        def find_out_base(self, y, x_B):
            # 找到x_B/y最小且y>0的位置
            index = []
            min_value = []
            for i, value in enumerate(y):
                if value <= 0:
                    continue
                else:
                    index.append(i)
                    min_value.append(x_B[i] / float(value))
    
            actual_index = index[np.argmin(min_value)]
            return actual_index
    
    
    if __name__ == "__main__":
        '''
        c = np.array([-20, -30, 0, 0]).reshape(-1, 1)
        A = np.array([[1, 1, 1, 0], [0.1, 0.2, 0, 1]])
        b = np.array([100, 14]).reshape(-1, 1)
        c = np.array([-4, -1, 0, 0, 0]).reshape(-1, 1)
        A = np.array([[-1, 2, 1, 0, 0], [2, 3, 0, 1, 0], [1, -1, 0, 0, 1]])
        b = np.array([4, 12, 3]).reshape(-1, 1)'''
        c = np.array([-3, -5, -4, 0, 0, 0]).reshape(-1, 1)
        A = np.array([[2, 3, 0, 1, 0, 0], [0, 2, 5, 0, 1, 0], [3, 2, 4, 0, 0, 1]])
        b = np.array([8, 10, 16]).reshape(-1, 1)
        simplex = Simplex(c, A, b)
        simplex.run()
    #转载自https://blog.csdn.net/cpluss/article/details/102596890

相关推荐