【机器学习5】python实现单纯形法和大M法

最新推荐文章于 2021-10-23 17:51:10 发布
最新推荐文章于 2021-10-23 17:51:10 发布
阅读量1.6k 收藏 21
点赞数 10
分类专栏： 人工智能
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/qq_34591921/article/details/105615957
版权

目录


写在前面：
关于单纯形法和大M法的原理，自己去百度吧，本篇文章只是给出了算法过程～

1.手算–单纯形法

题目：

在这里插入图片描述

手写计算最优解为为：27500

在这里插入图片描述

手算全过程请留言~

2.Python–单纯形法

题目依然是上面的↑

import numpy as np
def pivot(d,bn):
    l = list(d[0][:-2])
    jnum = l.index(max(l)) #转入编号
    m = []
    for i in range(bn):
        if d[i][jnum] == 0:
            m.append(0.)
        else:
            m.append(d[i][-1]/d[i][jnum])
    inum = m.index(min([x for x in m[1:] if x!=0]))  #转出下标
    s[inum-1] = jnum
    r = d[inum][jnum]
    d[inum] /= r
    for i in [x for x in range(bn) if x !=inum]:
        r = d[i][jnum]
        d[i] -= r * d[inum]        
def solve(d,bn):
    flag = True
    while flag:
        if max(list(d[0][:-1])) <= 0: #直至所有系数小于等于0
            flag = False
        else:
            pivot(d,bn)            
def printSol(d,cn):
    for i in range(cn - 1):
        if i in s:
            print("x"+str(i)+"=%.2f" % d[s.index(i)+1][-1])
        else:
            print("x"+str(i)+"=0.00")
    print("objective is %.2f"%(-d[0][-1]))
d = np.loadtxt("D:\\my.txt", dtype=np.float)
(bn,cn) = d.shape
s = list(range(cn-bn,cn-1)) #基变量列表
solve(d,bn)
printSol(d,cn)

计算结果：
在这里插入图片描述

3.Python包–单纯形法

题目依然是上面的↑

#导入包
from scipy import optimize
import numpy as np
#确定c,A_ub,B_ub
c = np.array([50,100])
A_ub = np.array([[1,1],[2,1],[0,1]])
B_ub = np.array([300,400,250])
#求解
res =optimize.linprog(-c,A_ub,B_ub)
print(res)

计算结果：
在这里插入图片描述

4.Excel–大M法

题目如下：

在这里插入图片描述

excel得到的结果如图：最优解为2

在这里插入图片描述

过程和上面的单纯法差不多 ~不清楚的先去学习单纯法吧 ~

5.python–大M法

题目是第4步骤贴出来的题目，代码参考: https://blog.csdn.net/desirew/article/details/79795554

导入txt文件内容格式为：
在这里插入图片描述

import os
import re
import numpy as np
import time
"""
输入线性规划模型到一个文档中，程序读取线性规划模型，化其为标准型。
找到初始可行基，建立单纯形表，迭代求解，判断解的类型，输出结果。
"""
#读取线性规划模型，并返回字符串列表
def inPut():
 
    #  读取线性规划模型
    model_str = []
    with open("大M法数据.txt",'r') as f:       #建立一个文本文档  储存线性规划模型
        lines = f.readlines()
    for element in lines:
        element = element[:-1]
        model_str.append(element)
    return(model_str)
 
#得到价值向量c
def find_C(model_str):
    x_temple =  re.findall(r"-?\d+\.?\d*",model_str[-1])
    x_temple = [int(i) for i in x_temple]
    x_max = max(x_temple)
    c_temple =  re.findall(r"-?\d+\.?\d*",model_str[0])
    c_temple = [float(i) for i in c_temple]
    c_list_1 = c_temple[1::2]
    c_list_2 = c_temple[::2]
    c=[]
    for i in range(x_max+1):
        c.append(0)
    for i in range(len(c_list_1)):
        c[int(c_list_1[i])] = c_list_2[i]
    c = np.array(c)
    return c
 
#求系数矩阵+资源向量
def find_A(model_str):
    b = []
    x_temple =  re.findall(r"-?\d+\.?\d*",model_str[-1])
    x_temple = [int(i) for i in x_temple]
    x_max = max(x_temple)
    A = [[0 for i in range(x_max+1)] for i in range(len(model_str)-2)]
    k = 0
    for element in model_str[1:-1]:
        x_temple =  re.findall(r"-?\d+\.?\d*",element)
        x_temple = [float(i) for i in x_temple]
        x_list_1 = x_temple[1::2]
        x_list_2 = x_temple[:-1:2]
        x=[]
        b.append(x_temple[-1])
        for i in range(x_max+1):
            x.append(0)
        for i in range(len(x_list_1)):
            x[int(x_list_1[i])] = x_list_2[i]
        for num in range(len(x)):
            A[k][num] = x[num]
        k += 1
    for num in range(len(b)):
        A[num][0] = b[num]
    A = np.array(A,dtype = np.float_)
    return A
 
#由以上三个函数可得到 1.系数+资源向量矩阵 2.价值向量    A c b  三个变量
 
#为进行资源向量非负化，约束条件等式化，记录每个约束的形式(  =记为0，>记为1,<记为-1)
def constraintCondition(model_str):
    yueshu = []
    for element in model_str[1:-1]:
        if '<' in element:
            yueshu.append(-1)
        elif '>' in element:
            yueshu.append(1)
        else:
            yueshu.append(0)
    yueshu = np.array(yueshu)
    return yueshu
 
#对A,b进行标准化
def change_A(str,A,c,yueshu):
    #约束条件等式化
    for i in range(A.shape[0]):
        if  A[i,0]< 0:
            yueshu[i]*=-1
            A[i]*=-1
    for i in range(A.shape[0]):
            if yueshu[i] == 1:
                #添加松弛变量
                temple_list = [0 for j in range(A.shape[0])]
                temple_list[i] = -1
                A = np.column_stack((A,temple_list))
 
    for i in range(A.shape[0]):
        if yueshu[i] == -1:
            #添加剩余变量
            temple_list = [0 for j in range(A.shape[0])]
            temple_list[i] = 1
            A = np.column_stack((A,temple_list))
        elif yueshu[i] == 1:
            #添加人工变量
            temple_list = [0 for j in range(A.shape[0])]
            temple_list[i] = 1
            A = np.column_stack((A,temple_list))
        elif yueshu[i] == 0:
            #  添加人工变量
            temple_list = [0 for j in range(A.shape[0])]
            temple_list[i] = 1
            A = np.column_stack((A,temple_list))
    return A
 
#对价值向量c进行标准化
def change_c(str,A,c,yueshu):
    if 'in' in str[0]:
        c = c*-1
    #约束条件等式化
    for i in range(A.shape[0]):
        if  A[i,0]< 0:
            yueshu[i]*=-1
 
    for i in range(A.shape[0]):
            if yueshu[i] == 1:
                #添加松弛变量
                c = list(c)
                c.append(0)
                c = np.array(c)
 
    for i in range(A.shape[0]):
        if yueshu[i] == -1:
            #添加剩余变量
            c = list(c)
            c.append(0)
            c = np.array(c)
 
        elif yueshu[i] == 1:
            #添加人工变量
            c = list(c)
            c.append(-pow(10,9))
            c = np.array(c)
        elif yueshu[i] == 0:
            #  添加人工变量
            c = list(c)
            c.append(-pow(10,9))
            c = np.array(c)
    return c
 
#基变量对应的价值系数
def C_B(x_b,c,m):
    c_b = []
    for element in x_b:
        c_b.append(c[element])
    return c_b
 
#计算检验数
def CN(A,c,x_b,c_b,m,n):
    cn = [0]
    for num_1 in range(1,n):
        cj = c[num_1]
        for num_2 in range(m):
            cj = cj - c_b[num_2] * A[num_2,num_1]
        if 0 < cj and cj < 0.0001:
            cj = 0
        cn.append(cj)
    return cn
 
#判断解的情况
def judge(A,cn,m,n,x_b,c_b):
    #最优解
    num = cn.count(0)
    if max(cn) <= 0:
        for j in range(len(c_b)):
            if c_b[j] == -pow(10,9):
                if A[j,0] !=0:
                    return 4                     #无可行解
        if num == m+1:
             return 1      #有最优解
        else:
            return 2      #有无穷多最优解
    else:
        for num in range(len(cn)):
            if cn[num] >0:
                lis = []
                for i in range(m):
                    lis.append(A[i,num])
                if max(lis) < 0:
                    for j in range(len(c_b)):
                        if c_b[j] == -pow(10,9):
                            if A[j,0] !=0:
                                return 4            #无可行解
                    return 3   #有无界解
    return 0      #继续迭代
 
#计算theta
def theTa(cn,A,m):
    theta = []
    j = cn.index(max(cn))
    for num in range(m):
        if A[num,j] > 0:
            theta.append(A[num,0] / A[num,j])
        else:
            theta.append(float("inf"))
 
    return theta
 
#进行基变换
def baseChange(A,cn,theta,m,n):
     # 找出中心元素
     j = cn.index(max(cn))
     i = theta.index(min(theta))
     main_elem = A[i,j]
 
     A[i] = A[i] / main_elem
     for num in range(m):
        if num != i:
            A[num] = A[num] - A[num][j] * A[i]
     return A
 
#进行基的换入与换出
def X_B(cn,theta,x_b):
    j = cn.index(max(cn))
    i = theta.index(min(theta))
    x_b[i] = j
    return x_b
 
#计算最优值
def Value(c_b,A,m):
    value = 0
    for num in range(m):
        value += c_b[num] * A[num,0]
    return value
 
#打印单纯形表
def Excel(A,c_b,x_b,c,theta,cn,m,n):
    print("%6s"%"Cj", end="")
    print("%6s"%"Cj", end="")
    print("%6s"%"Cj", end="")
    for i in c:
        if i != -pow(10,9):
            print("%6.1f"%i,end="")
        else:
            print("%6s"%"-M",end="")
    print("%6s"%"theta")
    #第二行
    print("%6s"%"CB", end="")
    print("%6s"%"XB", end="")
    print("%6s"%"b", end="")
    print("%6s"%"b", end="")
    for i in range(1,n):
        print("%6d" % i,end="")
    print("%6s"%"--")
 
    #打印数字
    for i in range(m):
        if c_b[i] != -pow(10,9):
            print("%6.1f" % c_b[i],end="")
        else:
            print("%6s" % "-M",end="")
        print("%6d" % x_b[i], end="")
        print("%6s" % "", end="")
        for j in range(n):
            print("%6.1f" % A[i,j] ,end="")
        print("%6.1f" % theta[i])
    #打印检验数
    print("%6s"%"CN", end="")
    print("%6s"%"CN", end="")
    print("%6s"%"CN", end="")
    for i in range(n):
        if cn[i]>10000 or cn[i]<-10000:
            print("%6s" %"MM", end="")
        else:
            print("%6.f" % cn[i],end="")
    print("%6s" % "--")
    print("")
 
 
def outPut(result,value,x_b,A,X):
    if result == 1:
        num = 0
        for i in x_b:
            if i <= len(X):
                if A[num,0]<0:
                    print("该问题无可行解:")
                    return
                X[i-1] = A[num,0]
            num +=1
        print("该问题有最优解:")
        print("X=(",end = "")
        for i in range(len(X)):
            if i == len(X)-1:
                print("%.1f"%X[i],end = "")
            else:
                print("%.1f,"%X[i],end = "")
        print(")")
        print("最优值:%.3f"%value)
    elif result == 2:
        print("该问题有无穷多解,其中一个最优解为:")
        num = 0
        for i in x_b:
            if i <= len(X):
                if A[num,0]<0:
                    return 4
                X[i-1] = A[num,0]
            num +=1
        print("X=(",end = "")
        for i in range(len(X)):
            if i == len(X)-1:
                print("%.0f"%X[i],end = "")
            else:
                print("%.0f,"%X[i],end = "")
        print(")")
        print("最优值:%.3f"%value)
 
    elif result ==3:
        print("该问题有无界解")
 
    elif result ==4:
        print("该问题无可行解")
 
 
 
#输入线性规划模型，找出变量
def solve():
    str = inPut()
    startTime = time.time()
    yueshu = constraintCondition(str)
    c = find_C(str)
    A = find_A(str)
    X = [ 0  for i in range(A.shape[1]-1)]
    #标准化
    A = change_A(str,A,c,yueshu)
    c = change_c(str,A,c,yueshu)
 
    #得到矩阵的维数
    m = A.shape[0]
    n = A.shape[1]
    #确定初始可行基X的下标,x_b存储可行基X的下标
    x_b = []
    for i in range(m):
        x_b.append(n-m+i)
 
    #X_B对应的价值系数c_b
    c_b = C_B(x_b,c,m)
    #初始单纯形表建立
 
    #1.计算检验数cn, theta
    #2.判断解的情况 若未达到终止条件 找到主元素 改变可行基 进行单位变换  进行第一步
    #3.若达到终止条件，输出结果
    num = 0
    while True:
        cn = CN(A,c,x_b,c_b,m,n)
        theta = theTa(cn,A,m)
        num += 1
        Excel(A,c_b,x_b,c,theta,cn,m,n)
        result = judge(A,cn,m,n,x_b,c_b)
        if result != 0:
            break
        theta = theTa(cn,A,m)
        A = baseChange(A,cn,theta,m,n)
        x_b = X_B(cn,theta,x_b)
        c_b = C_B(x_b,c,m)
    value = Value(c_b,A,m)
    if "in" in str[0]:
        value *= -1
    outPut(result,value,x_b,A,X)
    endTime = time.time()
    time_1 = endTime - startTime
    print("该问题迭代%d次,用时%.3fs" % (num,time_1))
 
 
solve()

运行结果：在这里插入图片描述
与excel计算结果一致！！

Ufoo360
关注
  • 10
    点赞
  • 踩
  • 21
    收藏
    觉得还不错? 一键收藏
  • 3
    评论
专栏目录
线性规划Python实现：使用库函数和不使用库函数进行单纯形法（大M法）线性规划
A981012的博客
04-21 1157
使用库函数和不使用库函数进行单纯形法线性规划使用scipy库进行单纯形法线性规划单纯形法大M法不使用库函数进行单纯形法线单纯形法大M法遇到的问题及解决 使用scipy库进行单纯形法线性规划 单纯形法 给出的线性规划为：max $2x_1 大M法 不使用库函数进行单纯形法线 单纯形法 大M法 遇到的问题及解决 ...
Python-单纯形法（大M法）求解 直接求解、借助scipy包
12-21
目录1、直接算法2、借助scipy库 在线性规划问题的约束条件中加人工变量后，要求在目标函数中相应地添加认为的M或一M为系数的项。在极大化问题中，对人工变量赋于一M作为其系数；在极小化问题中，对人工变量赋于一个M作为其系数，M为一任意大（而非无穷大）的正数。把M看作一个代数符号参与运算，是单纯形法求解的一种。 详细算法可参看小编的另一篇博客，Excel-单纯形法（大M法）求解 直接求解与规划求解功能。 下列代码的原题目如下所示： max 2×1+x2+x3 1、直接算法 代码如下所示： 如果没有安装numpy库的同学，请先 使用命令pip intall numpy安装。 # encoding
python实现单纯形法，大M法，拉格朗日乘子法
12-21
单纯形法： #导入包 from scipy import optimize import numpy as np #确定c,A,b,Aeq,beq c = np.array([115,90]) A = np.array([[10,20],[4,16],[15,10]]) b = np.array([200,128,220]) #Aeq = np.array([[1,-1,1]]) #beq = np.array([2]) #求解 res = optimize.linprog(-c,A,b) print(res) 输出结果： 大M法： #导入包 from scipy import opt
大M法的python编程求解和python包求解
qq_42585108的博客
04-19 1757
大M法的python编程求解和python包求解一、大M算法的求解步骤讲解二、python编程求解三、利用python包scipy的优化包optimize四、用excel求解五、分析结果 一、大M算法的求解步骤讲解 单纯形法的步骤是从一个初始极点出发，不断找到更优的相邻极点，直到找到最优的极点（或极线）。 消去xBxB x_BxB​得到问题的字典表达，即： mincTBB−1b+(cTN...
python实现单纯形法（大M法）
desirewen的博客
04-02 6534
import os import re import numpy as np import time """ 输入线性规划模型到一个文档中，程序读取线性规划模型，化其为标准型。 找到初始可行基，建立单纯形表，迭代求解，判断解的类型，输出结果。 """ #读取线性规划模型，并返回字符串列表 def inPut(): # 读取线性规划模型 model_str = [] wi...
Optimization:单纯形，最陡下降和黄金分割优化算法
03-06
在Python中，可以使用`scipy.optimize.linprog`库来实现单纯形法，解决线性目标函数在线性约束条件下的最小化问题。 其次，最陡下降法是梯度下降法的一种，用于求解无约束优化问题的全局最小值。该方法沿着目标函数...
机器学习 回归篇（1）——多元线性回归
12-20
机器学习 回归篇（1）——多元线性回归摘要线性回归简介python实现运行结果及可视化 摘要 本文介绍了最基础的回归问题——多元线性回归，并通过python进行实现及可视化展示运行结果。 线性回归简介 线性回归问题的...
simplex.py
04-03
最优化之线性规划求解算法之单纯性法的Python实现
通用单纯形——大M法程序
11-24
单纯形大M法求解最优化方法程序。简单适用！
单纯形法python
04-24
单纯形法python代码，单纯形法python代码，单纯形法python代码，单纯形法python代码，单纯形法python代码，单纯形法python代码，
最优化-线性规划单纯形法-大M法和两阶段法程序实现.rar
06-14
线性规划单纯形法-大M法和两阶段法程序实现（MATLAB）。程序有详细的注释。通俗易懂。 程序共有三个函数：main函数、twophase.m函数、MySimplex_method.m函数。 其中twophase.m函数是利用两阶段单纯形法。MySimplex_method.m函数是利用大M法单纯形法求解。 MySimplex_method.m函数既包含大M法，也包含单纯形法。twophase.m的单纯形法求解是通过调用MySimplex_method.m函数实现的。 在main函数中输入约束方程和目标函数，自动调用twophase.m函数、 MySimplex_method.m函数去进行规划求解。 如有需要，请联系私信我索取利用该程序求解实际问题的演讲PPT和Word。
线性规划大M法
01-28
大M法matlab代码 用于帮助学生学习在不能直接找到初始解的情况下如何使用大M法
Python库 | scikit_glpk-0.4.3-cp310-cp310-win32.whl
03-24
GLPK提供了各种算法，包括单纯形法、内点法等，能够处理大规模的线性模型，并且具有良好的可扩展性和灵活性。 **2. scikit-glpk与scikit-learn** scikit-glpk是针对Python的接口，使得Python开发者可以方便地利用...
Simplex单纯性算法的Python实现
modabobo
12-19 303
单纯性算法是解决线性规划的经典方法，上世纪50年代就提出了，其基本思想是在可行域内沿着边遍历所有的顶点，找出最优值，即为算法的最优值。 算法的执行过程如下： 求出初始基向量 构建单纯性表格 在所有非基向量对应的c中，找出一个最小的&lt;nobr style=""&gt;&lt;span class="math" id="MathJax-Span-71" style=""&gt;&amp
单纯形法中大m法_单纯形法（大M法）
weixin_28994403的博客
01-14 3102
function [x,fval,flag,iteration]=originalsimpleM(C,A,b)%原始单纯形法(大M法，无需给出初始基变量)%Programmed by Liyang(faruto's Studio~!) BNU MATH%Email:liyangbnu@mail.bnu.edu.cnQQ:516667408%last modified 2008.4.27%求解标准...
大M法的简单matlab程序
xuexijinhua的博客
10-23 2727
题目如上： 将其化为： min f=-50x1-40x1 -x1-x2<=-50 大M法函数代码如下： function [ A,b,c ] = M ( a,B,C) %使用前提：化为min z，Aixj之和<=0/=0，xj>0 M=10000; [m,n]=size(a); x=eye(m); A=[a,x]; b=B; c=[C,ones(1,m)*M]; %生成含M子块并加入子块到A和c中 i=1; [m,n]=size(A); for i=1:m ..
大M法
weixin_41443460的博客
04-26 3537
一、大M法的excel求解 1.1单纯形法 1.2大M法 1.2.1求解方法 1.2.2结果 求解列题： max 2x1+x2+x32x_1+x_2+x_32x1​+x2​+x3​ s.t={−2x2≥2x1−x2+x3=2x2−x3≤1x1,x2,x3≥0s.t=\begin{cases} -2x_2\geq 2 x_1-x_2+x_3= 2 x_2-x_3\leq 1 x_1,x...
机器学习梯度下降法python实现
最新发布
10-26
在Python中，可以使用NumPy和Matplotlib等库来实现梯度下降法。具体实现步骤如下： 1. 定义损失函数和梯度函数； 2. 初始化参数； 3. 设置学习率和迭代次数； 4. 循环迭代，更新参数，直到满足停止条件； 5. 可视化...

“相关推荐”对你有帮助么？

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 3
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值