如何解决ADMM 在凸 LP 上的实现,收敛问题
我在 cvxpy 上编写了一个简单的 LP,它有一个最佳答案。通过放宽一个约束,(p_p[0] + p_p[1] == 0),它可以作为两个独立的问题来解决,我已经使用 ADMM 解决了这个问题。这是在 for 循环中完成的,并行解决两个问题(prob_0 和 prob_1)并在每次迭代中更新对偶变量 (mu_bar)。出于某种我还没有弄清楚的原因,ADMM 实现在两个值上不断振荡,我没有得到收敛。您可以在下面的 colab 链接中看到代码: https://colab.research.google.com/drive/1WICekkYukrrezAVb1Yz-IP2hBfX3Qnpd?usp=sharing 任何有关如何解决该问题的见解都受到高度赞赏。谢谢。 ADMM 代码:
#%%%%%%%%%%%%%%%%%%%%%%% ADMM %%%%%%%%%%%%%%%%%%%%%
import cvxpy as cp
import numpy as np
import matplotlib.pyplot as plt
''' parameter and indices '''
p = 1
mp = 2
N = 12
NL = np.array([[-0.178625,-0.181325,-0.200575,-0.16895,-0.23155,-0.25675,-0.27425,-0.25575,-0.26,-0.312,-0.24,-0.24425 ],[3.007,2.875,2.884,2.839,3.104,2.853,2.889,2.771,2.803,2.718,3.018]])
###############################################################
'''decision variables'''
p_p = dict()
p_buy = dict()
p_sell = dict()
for i in range(mp):
p_p[i] = cp.Variable((N,p))
p_buy[i] = cp.Variable((N,p),nonneg = True)
p_sell[i] = cp.Variable((N,nonneg = True)
###############################################################
# dual variable,initiate with zero
mu_bar = dict()
mu_bar[0] = np.zeros((N,1))
# save primal variable of each problem
p0_val = dict()
p0_val[0] = np.zeros((N,1))
p1_val = dict()
p1_val[0] = np.zeros((N,1))
# save objective of each problem
obj0 = list()
obj1 = list()
# quadratic penalty coefficient
rho = 1
max_iter = 100
for i in range(max_iter):
# solve problem 0
obj_0 = cp.Minimize(
0.12 * 0.45 * sum(p_buy[0])
- 0.12 * 0.4005 * sum(p_sell[0])
+ 0.12 * 0.36 * sum(p_p[0])
+ mu_bar[i].T @ (p_p[0] + p1_val[i])
+ (rho/2) * cp.sum_squares(p_p[0] + p1_val[i])
)
constraints_0 = [ p_p[0] <= 10,p_p[0] >= -10,p_buy[0] <= 10,p_sell[0] <= 10,p_buy[0] - p_sell[0] + p_p[0] == np.array([NL[0,:]]).T,]
prob_0 = cp.Problem(obj_0,constraints_0)
prob_0.solve(solver = 'ECOS',verbose=False)
p0_val[i+1] = p_p[0].value
obj0.append(prob_0.value)
# solve problem 1
obj_1 = cp.Minimize(
0.12 * 0.45 * sum(p_buy[1])
- 0.12 * 0.4005 * sum(p_sell[1])
+ 0.12 * 0.36 * sum(p_p[1])
+ mu_bar[i].T @ (p_p[1] + p0_val[i])
+ (rho/2) * cp.sum_squares(p_p[1] + p0_val[i])
)
constraints_1 = [ p_p[1] <= 10,p_p[1] >= -10,p_buy[1] <= 10,p_sell[1] <= 10,p_buy[1] - p_sell[1] + p_p[1] == np.array([NL[1,]
prob_1 = cp.Problem(obj_1,constraints_1)
prob_1.solve(solver = 'ECOS',verbose=False)
p1_val[i+1] = p_p[1].value
obj1.append(prob_1.value)
p1_val[i+1] = np.array(p_p[1].value)
# update the dual variable
mu_bar[i+1] = mu_bar[i] + rho * (p0_val[i+1] + p1_val[i+1])
# plot the results
plt.figure(figsize=(14,4))
plt.title("check for (p_p[0] + p_p[1] = 0)")
plt.plot(p0_val[max_iter],label="p_p[0]")
plt.plot(p1_val[max_iter],label="p_p[1]")
plt.axhline(0,color='black',linestyle='--')
plt.legend()
plt.show()
plt.figure(figsize=(14,4))
plt.title("Objective functions")
plt.plot(obj0,label="OF_0")
plt.plot(obj1,label="OF_1")
plt.legend()
版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 dio@foxmail.com 举报,一经查实,本站将立刻删除。