|
楼主 |
发表于 2023-12-19 09:03:52
|
显示全部楼层
不使用scipy,帮我完善一下我的代码。给定A,b,c,能够输出无解、有无界解、有最优解的情况。
我的代码:
```
import numpy as np
def Simplex_eye(A, b, c):
"""
给定 max z =cx
s.t. Ax=b
x>=0
设A.shape=[m,n] (m<=n), r(A)= m,
b.shape=[m,1] 且b.all >= 0,
c.shape=[1,n],
x.shape=[n,1]
:param A: A中包含一个单位矩阵.
:param b:
:param c:
:return: x_opt:最优解。
fx_opt:最优值
iter:迭代次数
"""
m, n = A.shape
x_opt = []
fx_opt = []
iternum = 0
c = -1 * c # 把max问题变为min问题。
x = np.zeros(n)
# 1.找到A中的单位阵的列索引。
E_matrix = np.eye(m) # 创建单位矩阵
matching_columns = [] # 存储匹配的列数索引值
for i in range(m):
for j in range(n):
if np.array_equal(E_matrix[:, i], A[:, j]):
matching_columns.append(j)
break # E中的每个列向量能匹配到A中的一个就停下。防止匹配相同数值的列向量的索引。
# 迭代.
flag = True # 终止条件。
all_columns_idx = []
for i in range(n):
all_columns_idx.append(i)
while iternum < 10 and flag == True:
iternum += 1
remaining_columns = [i for i in all_columns_idx if i not in matching_columns] # 剩下的列索引是非基变量的索引。
# 2.最优性判别
B = A[:, matching_columns]
x[matching_columns] = np.dot(np.linalg.inv(B), b) # x_B
x[remaining_columns] = 0 # x_n
# 3.进基变量
f = np.dot(c[matching_columns], x[matching_columns])
w = np.dot(c[matching_columns], np.linalg.inv(B))
# 求判别数
dertas = np.dot(w, A[:, remaining_columns]) - c[remaining_columns]
if any(ele > 0 for ele in dertas):
x_in_idx = remaining_columns[dertas.argmax()] # 返回最大判别数对应的索引值。取索引i,得到进基变量xi。
else:
flag = False
x_opt = np.array(x)
fx_opt = f
break
# 4.离基变量
y = np.dot(np.linalg.inv(B), A[:, x_in_idx])
b_bar = np.dot(np.linalg.inv(B), b)
min_b_y = []
matching_out_idx = []
flag2 = True # 保证存在第一次的赋值。
for i, ele in enumerate(y):
if ele > 0:
b_y = (b_bar[i] / y[i])
if flag2:
min_b_y = b_y
matching_out_idx = i # 出基变量的索引
flag2 = False
if min_b_y > b_y:
min_b_y = b_y
matching_out_idx = i # 出基变量的索引
# if flag2 == True:
# iternum = 0
# break
# 5.换新基
matching_columns[matching_out_idx] = x_in_idx
return x_opt, fx_opt, iternum
if __name__ == "__main__":
# 定义系数矩阵 A 和约束向量 b, max优化问题的系数向量 c
A = np.array([[2, -3, 2, 1, 0],
[1/3, 1, 5, 0, 1]])
b = np.array([15, 20])
c = np.array([1, 2, 1, 0, 0])
# A = np.array([[-1, 2, 1, 0, 0],
# [ 2, 3, 0, 1, 0],
# [1, -1, 0, 0, 1]])
# b = np.array([4, 12, 3])
# c = np.array([4, 1, 0, 0, 0]) # 假如是max c,之后要在函数里面加负号。
# A = np.array([[1, -1, 1, 0], # 无界解
# [-3, 1, 0, 1]])
# b = np.array([2, 4])
# c = np.array([2, 3, 0, 0])
# A = np.array([[2, 7, 1, 0],
# [7, 2, 0, 1]])
# b = np.array([21, 21])
# c = np.array([4, 14, 0, 0])
# 调用。
x_opt, fx_opt, iternum = Simplex_eye(A, b, c)
x_opt = np.array([(round(x,4)) for x in x_opt])
fx_opt = round(fx_opt, 4)
# 输出结果
if iternum==10:
print('无解。')
else:
print("最优解x_opt = ", x_opt)
print("最优函数值fx_opt = ", fx_opt)
print("迭代次数为:", iternum)
``` |
|