Canethui 2020-03-03
MATLAB和Python解线性规划
//画可行域的方法 [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函数
使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);
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,[])
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