精华内容
下载资源
问答
  • 运筹系列1:线性规划单纯形法python代码

    万次阅读 多人点赞 2019-05-15 01:08:31
    单纯形法将搜索问题限制在了极点/极线范围内。那么如何找到它们呢?如有 xN=0x_N=0xN​=0,则 xxx就是极点了。那么此时有: 主问题可行条件: B−1b≥0B^{-1}b\ge 0B−1b≥0 主问题极点有很多,哪一个才是最优解呢...

    1. 线性规划问题

    线性规划标准模型如下:

    min ⁡ c x \min cx mincx
    s.t. A x = b Ax = b Ax=b
    x ≥ 0 x\ge 0 x0

    首先来看最简单的例子:A是一个方阵,比如:
    min ⁡ 2 x 1 + 5 x 2 \min 2x_1+5x_2 min2x1+5x2
    s.t. 3 x 1 + 3 x 2 = 6 3x_1+3x_2=6 3x1+3x2=6
    4 x 1 + x 2 = 5 4x_1+x_2=5 4x1+x2=5
    x 1 , x 2 ≥ 0 x_1,x_2 \ge 0 x1,x20

    只有一个确定解,最优值为7。

    约束条件不够的时候,表现为A变扁了。这时候的技巧就是,将A竖着切一刀,分离出一个方阵来:

    基变量:令 x = ( x B , x N ) x=(x_B,x_N) x=(xB,xN),其中 x B x_B xB满秩,称 x B x_B xB称为基变量 x N x_N xN称为非基变量

    将原问题用基变量重新表达为

    min ⁡ c B x B + c N x N \min c_Bx_B+c_Nx_N mincBxB+cNxN
    s.t. B x B + N x N = b Bx_B+Nx_N=b BxB+NxN=b
    x B ≥ 0 , x N ≥ 0 x_B\ge 0,x_N\ge 0 xB0,xN0

    极点:令 x x x m m m维(秩 m B ≤ m m_B\le m mBm),则需要 m m m个等式才能求解,它们的交点称为极点(有时候是极线)。

    单纯形法将搜索问题限制在了极点/极线范围内。那么如何找到它们呢?如有 x N = 0 x_N=0 xN=0,则 x x x就是极点了。那么此时有:

    主问题可行条件: B − 1 b ≥ 0 B^{-1}b\ge 0 B1b0

    主问题极点有很多,哪一个才是最优解呢?除了遍历搜索,还有一种方法:看对偶问题是否可行。上面问题的对偶问题为:

    max ⁡ b y \max by maxby
    s.t. y B T ≤ c B T yB^T\le c^T_B yBTcBT
    y N T ≤ c N T yN^T\le c^T_N yNTcNT

    B转置后仍然是满秩的方阵,因此等号是可以取到的,极点满足 y B T = c B T yB^T= c^T_B yBT=cBT,那么此时有:

    对偶问题可行条件: c B T B − 1 N ≤ c N T c^T_BB^{-1}N\le c^T_N cBTB1NcNT

    如果我们基向量拆分出来同时满足上面两个可行条件,那么恭喜找到最优解了:

    主问题可行:非基变量=0时,基变量系数 B − 1 b ≥ 0 B^{-1}b\ge 0 B1b0
    对偶问题可行:基变量系数=0时,非基变量系数 c N T − c B T B − 1 N ≥ 0 c^T_N-c^T_BB^{-1}N\ge 0 cNTcBTB1N0.
    定理:线性规划的极点同时满足主问题和对偶问题的可行条件时,这个极点是最优解。

    这两个条件非常重要!重要!重要!

    2. 求解步骤

    定理:线性规划问题的最优解一定是极点极线,并且满足比相邻极点都要优

    单纯形法的步骤是从一个初始极点出发,不断找到更优的相邻极点,直到找到最优的极点(或极线)

    消去 x B x_B xB得到问题的字典表达,即:
    min ⁡ c B T B − 1 b + ( c N T − c B T B − 1 N ) x N \min c^T_BB^{-1}b+(c^T_N-c^T_BB^{-1}N)x_N mincBTB1b+(cNTcBTB1N)xN
    s.t. B − 1 b − B − 1 N x N ≥ 0 B^{-1}b-B^{-1}Nx_N\ge 0 B1bB1NxN0
    x N ≥ 0 x_N\ge 0 xN0

    ( c N T − c B T B − 1 N ) i (c^T_N-c^T_BB^{-1}N)_i (cNTcBTB1N)i为第 i i i个非基变量的残差(reduced cost),如果残差小于0,那么说明这个非基变量还没有满足最优条件。直观上来看,这个非基变量取稍稍大于0的数就可以继续优化目标函数了。我们选择一个基变量和这个非基变量对换,就可以找到更优的相邻极点。
    另一种理解方法:残差对应的是对偶问题可行条件,小于0表示对偶问题还不可行,需要继续探索。

    接下来的问题是具体选择哪个基变量和哪个非基变量进行对换。我们用启发式原则,每次将负数残差 ( c N T − c B T B − 1 N ) i (c^T_N-c^T_BB^{-1}N)_i (cNTcBTB1N)i最小(绝对值最大)的非基变量 x i x_i xi替换为基变量,同时将 ( B − 1 b ( B − 1 N ) i ) j (\frac{B^{-1}b}{(B^{-1}N)_i})_j ((B1N)iB1b)j最小值对应的基变量 x j x_j xj替换为非基变量。这个进基/出基的过程称为pivoting
    另一种表达方式是: min ⁡ z = c x \min z = cx minz=cx,s.t. A x = b Ax = b Ax=b的pivoting是每次找出c中最小数对应的非基变量 x i x_i xi,再找出 b i / A i j b_i/A_{ij} bi/Aij最小的基变量 x j x_j xj进行对换。

    如果是max问题,令 c ′ = − c c'=-c c=c即可转化为min问题,相对应的,每次pivoting是找出 c ′ c' c最大值对应的非基变量。

    3. python算法实现

    这里假设原问题都是小于等于约束,这样添加松弛变量之后,问题一定有初始可行解;同时假设问题存在有限最优解。特殊情况将在下一节进行处理。代码为:

    import numpy as np   
    def solve():
        while max(d[-1][:-1]) > 0:
            jnum = np.argmax(d[-1][:-1]) #转入下标
            inum = np.argmin(d[:-1,-1]/d[:-1,jnum])  #转出下标
            s[inum] = jnum #更新基变量
            d[inum]/=d[inum][jnum]
            for i in range(bn):
                if i != inum:
                    d[i] -= d[i][jnum] * d[inum]
                
    def printSol():
        for i in range(cn - 1):
            print("x%d=%.2f" % (i,d[s.index(i)][-1] if i in s else 0))
        print("objective is %.2f"%(-d[-1][-1]))
    

    求解问题:max z = x 0 + 14 ∗ x 1 + 6 ∗ x 2 z = x_0+14*x_1+6*x_2 z=x0+14x1+6x2
    s.t.
    x 0 + x 1 + x 2 ≤ 4 x_0 + x_1 + x_2 \leq 4 x0+x1+x24
    x 0 ≤ 2 x_0 \leq 2 x02
    x 2 ≤ 3 x_2 \leq3 x23
    3 ∗ x 1 + x 2 ≤ 6 3*x_1 + x_2 \leq 6 3x1+x26

    添加变量将问题变为标准形,保存到data.txt中:

    1 1 1 1 0 0 0 4
    1 0 0 0 1 0 0 2
    0 0 1 0 0 1 0 3
    0 3 1 0 0 0 1 6
    1 14 6 0 0 0 0 0
    

    调用代码:

    d = np.loadtxt("data.txt", dtype=np.float)
    (bn,cn) = d.shape
    s = list(range(cn-bn,cn-1)) #基变量列表
    solve()
    printSol()
    

    运行后输出结果为:

    x0=0.00
    x1=1.00
    x2=3.00
    x3=0.00
    x4=2.00
    x5=0.00
    x6=0.00
    objective is 32.00
    

    4. 写后感

    将simplex用代码写出来,才觉得以前纠结那么久的问题原来那么简单。两三行代码能说清楚的事,何必写一堆看得人眼花缭乱的数学公式呢。
    另外,线性规划还有一些很基础的理论要掌握好:

    1. 极点和极方向的理论,这个是单纯型法的理论基础。可以参考这里
    2. 对偶理论,这个在以后经常会用到。
    展开全文
  • 单纯形法python

    2018-04-24 09:45:05
    单纯形法python代码单纯形法python代码单纯形法python代码单纯形法python代码单纯形法python代码单纯形法python代码
  • 单纯形法python实现

    千次阅读 2019-10-16 23:21:50
    单纯形法介绍 详见我的另一篇文章https://blog.csdn.net/cpluss/article/details/100806516 python代码 # coding=utf-8 # 单纯形法的实现,只支持最简单的实现方法 # 且我们假设约束矩阵A的最后m列是可逆的 # 这样就...

    单纯形法介绍

    详见我的另一篇文章https://blog.csdn.net/cpluss/article/details/100806516

    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)
        simplex = Simplex(c, A, b)
        simplex.run()
    
    

    只需要修改c、A、b便可以运行

    展开全文
  • 单纯形法讲解及Python代码实现一、了解单纯形法1.单纯形法的原理2.方法步骤二、例题讲解三、使用Python代码单纯形法求解线性规划最优解和最大值四、使用Python中scipy包进行上面的函数求解 一、了解单纯形法 1....
  • 本文为大家详细解读线性规划单纯形法的原理,以及超详细的python源码解读!

    1 单纯形法(Simplex method)

    1. 单纯形法是线性规划求解的经典方法之一,它的理论基础在于线性规划的解一定可以在顶点,即基可行解处取得。
    2. 举个简单例子,大家熟悉的二元线性规划中,目标函数的直线在坐标系上移动,在与可行域的交界处会取得最优解,即使目标函数与边界线相重合,也必然会有顶点位于该重合的直线上,即在顶点处必然会取得最优解。
    3. 基于以上论述,单纯形法就是从一个基可行解开始,通过一次改变一个基变量,从而实现从当前顶点转至下一个顶点的过程,结合基变量的模型表达式如下:
      m a x c B X B + c N X N max\quad c_B X_B+c_N X_N maxcBXB+cNXN
      s . t . B X B + N X N = b s.t.\quad B{X}_{B}+N{X}_{N}=b s.t.BXB+NXN=b
      X B ≥ 0 , X N ≥ 0 \qquad X_{B}\ge0,X_{N}\ge0 XB0,XN0

    2 编程思路

    对于单纯形法的编程实现,其主要基于单纯形表,即约束系数矩阵,约束常数项以及检验数组成,通过初等行变换不断构建列线性无关的基,观察检验数的变化,当检验数表示其他变量入基不会带来更好的目标时,也就达到了最优解,具体的过程解读以及内在含义将在下方结合代码说明。

    3 python实现原理解读

    鉴于博主也为算法小白,源代码参考于运筹系列1:线性规划单纯形法python代码
    原博主的代码实现方式非常简洁,但对于初学者而言还是会相对难以理解,本文也旨在对其做较为详细的说明与解读。

    1. 数据准备
      该算法的输入为一个已有一组基的单纯形表,目标函数取最大值,如下表

    ( 1 14 6 0 0 0 0 0 1 1 1 1 0 0 0 4 1 0 0 0 1 0 0 2 0 0 1 0 0 1 0 3 0 3 1 0 0 0 1 26 ( x 1 x 2 x 3 x 4 x 5 x 6 x 7 b ) ) \quad\quad\quad\quad\quad\quad\quad\quad\begin{pmatrix} 1&14&6&0&0&0 &0&0\\\\ 1&1&1&1&0&0 &0&4\\\\ 1&0&0&0&1&0 &0&2\\\\ 0&0&1&0&0&1 &0&3\\\\ 0&3&1&0&0&0 &1&26\\\\ (x_1&x_2&x_3&x_4&x_5&x_6 &x_7&b)\\\\ \end{pmatrix} 11100(x1141003x261011x301000x400100x500010x600001x7042326b)

    在该矩阵中,第一行(row)的含义为目标函数值由非基变量表达的式子,即 Z = σ j X j + Z 0 Z=\sigma_jX_j+Z_0 Z=σjXj+Z0,其第一行即为 Z = X 1 + 14 X 2 + 6 X 3 + 0 Z=X_1+14X_2+6X_3+0 Z=X1+14X2+6X3+0 注意,第一行最后一列为 Z 0 Z_0 Z0,由于该算例的初始基可行解均为辅助变量,即 X 4 , X 5 , X 6 , X 7 X_4,X_5,X_6,X_7 X4X5X6X7 ,则刚好原问题中的初始变量的 c j c_j cj都得以保留。(在一般问题当中也较为普遍,若每一行都添加了辅助变量来构造初始基可行解)

    由于此时目标函数由一系列非基变量的线性组合构成,即与基变量的取值无关,则可通过该表达式观察是否有更好的解,在此例中,可观察 X 2 X_2 X2的系数最大,即 X 2 X_2 X2增加一个单位,目标函数值即可增加14个单位,则 X 2 X_2 X2应当入基从而获得更优基。 没错,这一行的实际意义就是检验数 σ j = c j − ∑ i = 1 m c m a i j \sigma_j =c_j-\sum_{i=1}^mc_ma_{ij} σj=cji=1mcmaij, 即该非基变量增加一个单位,其收益 c j c_j cj与在约束条件下,基变量随其增加1而相应减少的收益 ∑ i = 1 m c m a i j \sum_{i=1}^mc_ma_{ij} i=1mcmaij之差,这个式子也好理解,比如看矩阵第五行,即第四个约束条件,若 X 2 X_2 X2增加一个单位,则为了满足约束,基变量 X 7 X_7 X7需减少3个单位,其收益损失 c 7 ∗ a 57 = 0 ∗ 3 = 0 c_7*a_{57}=0*3=0 c7a57=03=0个单位,其余基变量也同理,这也就是检验数的实际效果。

    根据以上推断,在编写代码时,我们要做到不断的变换基变量,直至:

    1. 第一行表达式中的非基变量如果入基,则其需要离开表达式,保证表达式中没有基变量,否则变换非基变量,基变量同样也会改变,则目标函数值变化方向不定,检验数失去效果。
    2. 第一行的前七个数均小于零,即满足最优性检验条件

    有了这样的准则,我们只需编写数据读取与初始化操作,出基入基操作(pivot),与最优性判断条件(嵌在solve方法中)。

    4 python代码

    1. 数据读取与初始化
      数据如下:
    
    1 14 6 0 0 0 0 0
    
    1 1 1 1 0 0 0 4
    
    1 0 0 0 1 0 0 2
    
    0 0 1 0 0 1 0 3
    
    0 3 1 0 0 0 1 6
    
    ##matrix即为上面展示的矩阵
    matrix = np.loadtxt("data.txt", dtype=float)
    (b_num, c_num) = matrix.shape
    # b_num为约束行数加第一行的目标函数表达式行,c_num为约束系数列加上常数项b列
    indexes_B = list(range(c_num - b_num, c_num - 1)) #基变量列表
    solve()
    printSol()
    
    
    1. 出基入基操作Pivot()
    def pivot():
        # 求转入基变量的列索引id
        l = list(matrix[0][:-1])
        index_intoB = l.index(max(l))#即取检验数的最大值列索引
    
        # 算出基变量在表中的行索引,用b列除以入基变量列后的值进行比较判断
        col_inB = []#入基的基变量的列,看哪个基变量出基
        for i in range(b_num):
            if matrix[i][index_intoB] == 0:
                col_inB.append(0.)
            else:
                col_inB.append(matrix[i][-1] / matrix[i][index_intoB])
        index_outB = col_inB.index(min([x for x in col_inB[1:] if x > 0]))#取b/aij的最小值,否则其他约束会不可行
    
        #将基变量列表中进行新旧基的替换
        indexes_B[index_outB - 1] = index_intoB #index_outB - 1为实际出基变量在基变量列表中的索引,因为矩阵第一行是c,约束矩阵A从第二行开始
        cur_cell_to1 = matrix[index_outB][index_intoB] #新基的应为1的位置当前的数值
        matrix[index_outB] /= cur_cell_to1 #该行全部除以该值,使其为1,即该行标准化
    
        #通过行变换将其他行非零元去除
        #注意,这一操作对第一行也执行了,是在将表达式中的新基变量也给消除!用其他非基变量来表示!第一行的最后一个数值极为z0,即目标函数值!
        for row in [x for x in range(b_num) if x != index_outB]:
            cur_cell_to0 = matrix[row][index_intoB]
            matrix[row] -= cur_cell_to0 * matrix[index_outB]
            
    
    1. 算法执行与最优性检验
    def solve():
        flag = True
        while flag:
            if max(list(matrix[0][:-1])) <= 0:  # 直至所有系数小于等于0
                flag = False
            else:
                pivot()
    
    def printSol():
        for i in range(c_num - 1):
            if i in indexes_B:
                print("x" + str(i) + "=%.2f" % matrix[indexes_B.index(i) + 1][-1])
            else:
                print("x" + str(i) + "=0.00")
        print("objective is %.2f" % (-matrix[0][-1]))
        
    
    1. 整体代码(考虑编译先后顺序)
    import numpy as np
    
    def pivot():
        # 求转入基变量的列索引id
        l = list(matrix[0][:-1])
        index_intoB = l.index(max(l))#即取检验数的最大值列索引
    
        # 算出基变量在表中的行索引,用b列除以入基变量列后的值进行比较判断
        col_inB = []#入基的基变量的列,看哪个基变量出基
        for i in range(b_num):
            if matrix[i][index_intoB] == 0:
                col_inB.append(0.)
            else:
                col_inB.append(matrix[i][-1] / matrix[i][index_intoB])
        index_outB = col_inB.index(min([x for x in col_inB[1:] if x > 0]))#取b/aij的最小值,否则其他约束会不可行
    
        #将基变量列表中进行新旧基的替换
        indexes_B[index_outB - 1] = index_intoB #index_outB - 1为实际出基变量在基变量列表中的索引,因为矩阵第一行是c,约束矩阵A从第二行开始
        cur_cell_to1 = matrix[index_outB][index_intoB] #新基的应为1的位置当前的数值
        matrix[index_outB] /= cur_cell_to1 #该行全部除以该值,使其为1,即该行标准化
    
        #通过行变换将其他行非零元去除
        #注意,这一操作对第一行也执行了,是在将表达式中的新基变量也给消除!用其他非基变量来表示!第一行的最后一个数值极为z0,即目标函数值!
        for row in [x for x in range(b_num) if x != index_outB]:
            cur_cell_to0 = matrix[row][index_intoB]
            matrix[row] -= cur_cell_to0 * matrix[index_outB]
    
    def solve():
        flag = True
        while flag:
            if max(list(matrix[0][:-1])) <= 0:  # 直至所有系数小于等于0
                flag = False
            else:
                pivot()
    
    def printSol():
        for i in range(c_num - 1):
            if i in indexes_B:
                print("x" + str(i) + "=%.2f" % matrix[indexes_B.index(i) + 1][-1])
            else:
                print("x" + str(i) + "=0.00")
        print("objective is %.2f" % (-matrix[0][-1]))
    
    matrix = np.loadtxt("data.txt", dtype=float)
    (b_num, c_num) = matrix.shape
    # b_num为约束行数加第一行的目标函数表达式行,c_num为约束系数列加上常数项b列
    indexes_B = list(range(c_num - b_num, c_num - 1)) #基变量列表
    solve()
    printSol()
    
    

    5 后记

    博主研究方向是高速铁路运行图(大规模混合整数规划)的优化,现在正在啃整数规划的基础原理,这个单纯形法是我研学之路的开始,未来会对各类整数规划经典问题进行学习并分享我的经验,欢迎大家交流!

    展开全文
  •   本文所提供的单纯形法Python实现,基于sympy和numpy库。使用时请先安装相关库。   优点:可以直接按目标函数和不等式约束的原形式输入。   缺点(BUG): 所有变量必须>=0 未解决约束条件全为等式情况   ...
  • python实现线性规划中的单纯形法 例题如下(已是标准形式): maxz=1500x1+1000x2max z=1500x_1+1000x_2maxz=1500x1​+1000x2​ {3x1+2x2+x3=652x1+x2+x4=403x2+x5=75x1,x2,x3,x4,x5≥0 \left\{ \begin{array}{c} ...

    用python实现线性规划中的单纯形法

    例题如下(已是标准形式):
    m a x z = 1500 x 1 + 1000 x 2 max z=1500x_1+1000x_2 maxz=1500x1+1000x2

    { 3 x 1 + 2 x 2 + x 3 = 65 2 x 1 + x 2 + x 4 = 40 3 x 2 + x 5 = 75 x 1 , x 2 , x 3 , x 4 , x 5 ≥ 0 \left\{ \begin{array}{c} 3x_1+2x_2+x_3=65 \\ 2x_1+x_2+x_4=40 \\ 3x_2+x_5=75\\ x_1,x_2,x_3,x_4,x_5\geq0 \end{array} \right. 3x1+2x2+x3=652x1+x2+x4=403x2+x5=75x1,x2,x3,x4,x50

    代码如下:

    #输入变量的数量
    arg_num = int(input())
    #输入目标函数
    max_z = list(map(int,input().split(' ')))
    z = 0
    x = {}
    
    class pure_chat():
        def __init__(self):
            self.c_b = []
            self.b = []
            self.opt = []
            self.no_base_args = []
            self.ex_num_list = max_z[:]
    	
    	#输入约束条件
        def input_chat(self):
            for i in range(len(max_z)):
                if max_z[i] != 0:
                    self.no_base_args.append(i)
    
            while True:
                exper = list(map(int, input().split(' ')))
                if exper.count(0) == arg_num + 1:
                    break
                self.opt.append(exper)
    
            for i in range(arg_num - len(self.opt) + 1):
                self.b.append(arg_num - len(self.opt) + i)
                self.c_b.append(max_z[arg_num - len(self.opt) + i])
    
        def get_max(self):
            m_i = 0
            max = -10000
            for i in self.no_base_args:
                if self.ex_num_list[i] > max:
                    max = self.ex_num_list[i]
                    m_i = i
            return m_i
    
    	#判断非基变量的检验数是否都为小于0
        def examine(self):
            for i in self.no_base_args:
                if self.ex_num_list[i] > 0:
                    return False
            return True
            
    	#确定出基变量
        def out_base(self):
            min = 100000
            i = self.get_max()
            if self.ex_num_list[i] > 0:
                for j in range(len(self.opt)):
                    if self.opt[j][i] == 0:
                        continue
                    t = self.opt[j][arg_num] / self.opt[j][i]
                    if t < min:
                        min = t
                        out_i = i
                        out_j = j
            return out_i, out_j
    	
    	#将基变量矩阵化为单位矩阵
        def to_e(self,b_i, b_j):
            s = self.opt[b_i][b_j]
            for i in range(len(self.opt[b_i])):
                self.opt[b_i][i] = self.opt[b_i][i] / s
            for j in range(len(self.opt)):
                if j != b_i:
                    t = self.opt[j][b_j]
                    for k in range(len(self.opt[j])):
                        self.opt[j][k] = self.opt[j][k] - t * self.opt[b_i][k]
    
    	#出基与进基
        def change_base_args(self,b_i, b_j):
            t = self.no_base_args.index(b_j)
            self.no_base_args[t] = self.b[b_i]
            self.b[b_i] = b_j
    
    	#更新基变量的系数矩阵
        def update_c_b(self):
            for i in range(len(self.b)):
                self.c_b[i] = max_z[self.b[i]]
    	
    	#计算检验数
        def compute_ex_num(self,b_i, b_j):
            self.ex_num_list[b_j] = 0
            for i in self.no_base_args:
                s = max_z[i]
                count = 0
                for j in self.b:
                    s = s - max_z[j] * self.opt[count][i]
                    count += 1
                self.ex_num_list[i] = s
    
    	#计算目标函数值
        def compute_z(self):
            s = 0
            for i in range(len(self.opt)):
                s = s + self.c_b[i] * self.opt[i][arg_num]
            return s
    
    	#计算最优解
        def compute_x(self):
            for i in range(len(self.opt)):
                x[self.b[i]] = self.opt[i][arg_num]
            for i in range(arg_num):
                if i not in self.b:
                    x[i] = 0
    
    
    chat = pure_chat()
    chat.input_chat()
    
    #循环迭代
    while True:
        if chat.examine():
            chat.compute_x()
            x = sorted(x.items(), key=lambda x: x[0])
            z = chat.compute_z()
            break
        j,i = chat.out_base()
        chat.change_base_args(i,j)
        chat.update_c_b()
        chat.to_e(i,j)
        chat.compute_ex_num(i,j)
    
    print(x)
    print(z)
    
    '''
    
    测试数据:
    5
    1500 1000 0 0 0
    3 2 1 0 0 65
    2 1 0 1 0 40
    0 3 0 0 1 75
    0 0 0 0 0 0
    
    输出的结果:
    [(0, 15.0), (1, 10.0), (2, 0), (3, 0), (4, 45.0)]
    32500.0
    
    '''
    
    

    输入的数据格式为:
    5
    1500 1000 0 0 0
    3 2 1 0 0 65
    2 1 0 1 0 40
    0 3 0 0 1 75
    0 0 0 0 0 0

    第一行的5为变量的个数,第二行为目标函数中的对应系数,其余行为约束条件的对应系数,最后输入全0表示结束输入

    输出的结果为:
    [(0, 15.0), (1, 10.0), (2, 0), (3, 0), (4, 45.0)]
    32500.0

    第一行的列表表示最优解,为 x 1 = 15 , x 2 = 10 , x 3 = 0 , x 4 = 0 , x 5 = 45 x_1=15,x_2=10,x_3=0,x_4=0,x_5=45 x1=15,x2=10,x3=0,x4=0,x5=45;第二行为最优解对应的最优目标函数值,为32500.

    ps:个人github地址:https://github.com/zhouyanb
    github地址

    展开全文
  • 1,python 代码 这里直接附上代码: #!/usr/bin/python # coding:utf-8 # 19-3-20 下午3:12 # @File : Simplex.py import numpy as np # import os data= [] #全局变量用来存储矩阵 def pivot(...
  • 单纯形法求解线性规划问题python代码 输入格式代码中可见 #simplex # #单纯形法求解一般形式的线性规划 # x1+x2<5 # 3x1+x2<8 # x1,x2>=0 # max(2x1+x2) #初始单纯形表 # x1 x2 x3 x4 #x3 1 1 1 0 5 #x4 3...
  • 单纯形法——Python实现(一)

    千次阅读 2020-04-03 16:38:42
    simplex.py 代码内容参见 我的上传 test.py 代码内容如下: from sympy import * from simplex import * def example_0(): x1, x2, x3 = symbols('x1, x2, x3') obj = -3 * x1 + x2 + x3 variables = [x1, x2...
  • Python求解线性规划问题_两阶段法实现的单纯形法,包括.py和.ipynb两种格式,用Jupyter Notebook打开.ipynb或者用Python软件打开.py都可成功运行,压缩包中包括测试数据,代码可输出唯一解,无穷多解,无界解,无解...
  • 人工智能-线性规划(单纯形法、大M法)和非线性规划(拉格朗日乘子法) 一、实验内容: 二、相关算法介绍 1、线性规划 线性规划(Linear programming,简称LP),是运筹学中研究较早、发展较快、应用广泛、方法较...
  • Python实现单纯形法

    2021-01-04 15:20:34
    import tkinter import tkinter.filedialog import tkinter.messagebox import tkinter.ttk ...root.title('单纯形法') root['height']=250 root['width']=420 #放'目标函数:' labelName=tkinter.Label(root,
  • 优化方法总结续篇:下降单纯形法(downhill simplex) 及python示例代码 下降单纯形法(downhill simplex method)是一个广泛使用的“derivative free”的优化算法。一般来说它的效率不高,但是文献[1]...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 1,230
精华内容 492
关键字:

单纯形法python代码

python 订阅