#利用Python进行矩阵运算
impor numpy as np

matrix_1 = np.array([[1, 2, 3], [2, 3, 4], [3, 4, 5], [4, 5, 6]])   #这是一个四行三列的矩阵，其列表是竖着放置的
matrix_2 = np.array([[i for i in range(3)] for x in range(4)])
matrix_3 = np.array(np.arange(12).reshape([4, 3]))
print(matrix_1 + matrix_2)  #两个矩阵相加（对应位置相加）
print(matrix_1 * matrix_2)  #矩阵的元素乘法（对应位置相乘）
print(matrix_3.dot(matrix_1.T)) #元素的矩阵乘法（n行乘以n列加到第n行n列中）；矩阵.T表示矩阵的转置
matrix = np.array([[1, 2], [2, 1]])
result = np.linalg.inv(matrix)  #求一个矩阵的逆矩阵（前提这个矩阵可逆，不可逆会报错）
print(result)
print(np.linalg.det(matrix))    #求一个矩阵的行列式（最后的结果是一个数字）
print(np.linalg.matrix_rank(matrix))    #求一个矩阵的秩

#使用Python求解线性函数
'''
10x-y-2z=72
-x+10y-2z=83
-x-y+5z=42
'''
import numpy as np

coefficient_matrix = np.array([[10, -1, -2], [-1, 10, -2], [-1, -1, 5]])    #这是线性方程的系数矩阵
result = np.array([72, 83, 42]) #这是其右边的常数
inv_matrix = np.linalg.inv(coefficient_matrix)  #系数矩阵的逆
x = inv_matrix.dot(result)    #系数矩阵的逆与常数做点积，得出有关解的矩阵
print(x)
x = np.linalg.solve(coefficient_matrix, result) #对于上面的26,27行，可以用29行代替，返回的还是解的列表
print(x)

    #同样可以用sympy库来实现线性计算

from sympy import symbols, Eq, solve

x, y, z = symbols('x, y, z')
eqs = [
    Eq(10*x-y-2*z, 72),
    Eq(-x+10*y-2*z, 83),
    Eq(-x-y+5*z, 42)
]
result = solve(eqs, [x, y, z])
print(result)

求解线性规划问题

'''线性规划的标准形式
    z = min cx
    A_ub . x <= b_ub
    Aeq.x = Beq
    l <= x <= u bouds[l, u]没有则填null
'''

'''求下列线性规划
    z = max 2x + 3y - 5z    #求最大值
    x + y + z = 7   #下面三个都是限制条件
    2x - 5y + z >= 10
    x + 3y + z <= 12
    x, y, z >= 0
'''

#标准形式是求最小值，要求最大值就将矩阵都加符号，标准形式是小于等于，限制条件是大于等于就将矩阵加个符号
c = np.array([-2, -3, 5])  #由于要求的是最大值，因此将矩阵改为负号
Aeq = np.array([[1, 1, 1]]) #Aeq和A_ub一定要是一个二维的
Beq = np.array([7])
A_ub = np.array([[-2, 5, -1], [1, 3, 1]])
B_ub = np.array([-10, 12])
x = [0, None]
y = (0, None)
z = [0, None]
result = scipy.optimize.linprog(c, A_ub, B_ub, Aeq, Beq, [x, y, z])
print(result)   #输出中fun为要求的最值，x为最优解

#使用PulP来实现线性规划求解

import pulp
#首先创立一个线性规划问题
my_LPproblem = pulp.LpProblem("my_problem", sense=pulp.LpMaximize) #第一个参数为问题的名称（可以不写），sense关键字参数为求解最大值还是最小值
#设置变量，第一个参数是变量名称，后面分别是参数的范围
x = pulp.LpVariable("x", 0, None)
y = pulp.LpVariable("y", 0, None)
z = pulp.LpVariable("z", 0, None)
#加入最值函数(objective属性)
my_LPproblem += 2*x + 3*y - 5*z
#加入约束方程
my_LPproblem += x + y + z == 7
my_LPproblem += 2*x - 5*y + z >= 10
my_LPproblem += x + 3*y + z <= 12
#求解
result = my_LPproblem.solve()
print(result)
#打印取得最值后的变量
print(f"x={x.varValue}, y={y.varValue}, z={z.varValue}")
#打印最值
print(pulp.value(my_LPproblem.objective))

解决运输类问题

#使用Python解决运输问题
import numpy as np
import pulp

def Transportation_problem(cost, x_max, y_max):
    y_len = len(cost)   #求解x，y轴的长度
    x_len = len(cost[0])
    #建立问题
    my_problem = pulp.LpProblem("Transportation_problem", sense=pulp.LpMaximize)

    #这个是用于创建变量Xi,j，cat传参表示取整数
    var = [[pulp.LpVariable(f"x{y}{x}", lowBound=0, upBound=None, cat=pulp.LpInteger) for x in range(x_len)] for y in range(y_len)]

    #这个函数的作用相当于将一个多维数组扁平化，变成一维数组，这里如果使用的是numpy中的ndarray对象可以直接调用flatten方法(使用下面两条语句也可以达到同样的效果)
    #ndvar = np.array(var)
    #my_problem += pulp.lpDot(ndvar.flatten(), cost.flatten())
    flatten = lambda x: [y for l in x for y in flatten(l)] if type(x) is list else [x]
    my_problem += pulp.lpDot(flatten(var), cost.flatten())  #将两个一维数组相乘，得出最值函数的表达式(价格函数)，加入到问题中

    #每一行相加应该小于这一行的总和的最大值
    for y in range(y_len):  #每一行上的约束条件
        my_problem += pulp.lpSum(var[y]) <= x_max[y]
    #每一列相加小于这一列的最大值
    for x in range(x_len):  #每一列上的约束条件
        my_problem += pulp.lpSum([var[y][x] for y in range(y_len)]) <= y_max[x]

    my_problem.solve()
    #返回最优解和取得最优解时的分配
    return {"objective": pulp.value(my_problem.objective), "var": [[pulp.value(var[y][x]) for y in range(y_len)] for x in range(x_len)]}

cost = np.array(
    [
        [500, 550, 630, 1000, 800, 700],
        [800, 700, 600, 950, 900, 930],
        [1000, 960, 840, 650, 600, 700],
        [1200, 1040, 980, 860, 880, 780],
    ]
)
x_max = np.array([76, 88, 96, 40])
y_max = np.array([42, 56, 44, 39, 60, 59])
result = Transportation_problem(cost, x_max, y_max)
print(result)

整数约束问题

import numpy as np
import math
import scipy
import sys

def IntegerProblem(c, A_ub, B_ub, A_eq, B_eq):
    tmin = 1.0E-12  #设置爱一个很小的数，为了后面判断整数使用

    #先求解这个线性规划
    result = scipy.optimize.linprog(c=c, A_ub=A_ub, b_ub=B_ub, A_eq=A_eq, b_eq=B_eq)

    bestX = [sys.maxsize]*len(c)    #sys.maxsize是一个很大的值，这里创建一个和未知数一样大的数组
    if not (type(result.x) is float or result.status != 0):
        bestX = result.x    #如果不是浮点数就可以直接读取答案

    #这个是最优解
    bestVal = sum([x*y for x, y in zip(c, bestX)])

    #floor和ceil方法分别是取比一个小数小的最大整数和比一个小数大的最小整数，这里要想不小于，只有x为一个整数（浮点型整数？）
    if all(((x - math.floor(x) < tmin) or math.ceil(x) - x < tmin) for x in bestX):
        return [bestX, bestVal]
    else:   #分支定界
        ind = [i for i, x in enumerate(bestX) if x - math.floor(x) > tmin and math.ceil(x) - x > tmin][0]   #后面两个条件用来判断哪个位置的x是浮点数，enumerate返回的是某个元素的位置和其值
        newCon1 = [0]*len(A_ub[0])  #相当于创建了线性规划的不等式的右边，但是参数都是0
        newCon2 = [0]*len(A_ub[0]) #同样，这两步是为了后面添加条件x>=floor(x)和x<=ceil(x)，这两个后面放在A_ub中判断到底是大好还是小号
        newCon1[ind] = -1   #这个是大于号，所以左右两边都取负号(x>=floor(x))
        newCon2[ind] = 1    #条件x<=ceil(x)
        new_A1_ub = A_ub.copy()
        new_A2_ub = A_ub.copy()
        new_A1_ub.append(newCon1)   #添加限制条件x>=ceil(x)
        new_A2_ub.append(newCon2)   #另一限制条件x<=floor(x)，比较两者的大小
        new_B1_ub = B_ub.copy()
        new_B2_ub = B_ub.copy()
        new_B1_ub.append(-math.ceil(bestX[ind]))   #添加ceil(x)，新的x的取值大于等于比x大的最小整数
        new_B2_ub.append(math.floor(bestX[ind]))

        #下面使用递归重新计算了不同的添加的两种限制条件下的极值，返回优的那一个
        result_1 = IntegerProblem(c, A_ub=new_A1_ub, B_ub=new_B1_ub, A_eq=A_eq, B_eq=B_eq)
        result_2 = IntegerProblem(c, A_ub=new_A2_ub, B_ub=new_B2_ub, A_eq=A_eq, B_eq=B_eq)

        if result_1[1] > result_2[1]:   #比较两者的bestVal
            return result_2
        else:
            return result_1

c = [3, 4, 1]
A = [[-1, -6, -2], [-2, 0, 0]]
b = [-5, -3]
Aeq = [[0, 0, 0]]
beq = [0]
print(IntegerProblem(c, A, b, Aeq, beq))

非线性规划问题

#非线性规划

    #例子：计算1/x+x的最小值
import numpy as np
import scipy

def fun(args):
    a = args
    v = lambda x: a/x[0]+x[0]   #这里的x[]后面的标是指第几个未知数，比如第一个未知数就是x[0]第二个未知数就是x[1]
    return v

args = (1)  #这里应该传递元组
x0 = np.asarray((2))    #这里是猜想值
result = scipy.optimize.minimize(fun(args), x0, method="SLSQP") #这个得出来的只是局部最优解
print(result)

#计算(2+x[1])/(1+x[2])-3x[1]+4x[3]，其中x[i]的范围都在0.1-0.9中间

import scipy

def fun(args):
    a, b, c, d = args
    v = lambda x:(a + x[0]) / (b + x[1]) - c*x[0] + 4*x[2]
    return v

def con(args):  #约束条件，分为eq和ineq，eq表示后面算式结果等于0，ineq表示后面算式结果大于0
    x0min, x0max, x1min, x1max, x2min, x2max = args
    cons = (
        {"type":"ineq", "fun":lambda x: x[0] - x0min},
        {"type":"ineq", "fun":lambda x: x0max - x[0]},
        {"type":"ineq", "fun":lambda x: x[1] - x1min},
        {"type":"ineq", "fun":lambda x: x1max - x[1]},
        {"type":"ineq", "fun":lambda x: x[2] - x2min},
        {"type":"ineq", "fun":lambda x: x2max - x[2]}
    )
    return cons

args = (2, 1, 3, 4)
ranges = (0.1, 0.9, 0.1, 0.9, 0.1, 0.9)
cons = con(ranges)
x = np.asarray((0.5, 0.5, 0.5))
result = scipy.optimize.minimize(fun(args), x, method="SLSQP", constraints=cons)
print(result.fun)   #最优解
print(result.success)   #是否达到容忍值范围
print(result.x) #局部最优解时参数的取值