微信公众号搜"智元新知"关注
微信扫一扫可直接关注哦!

ADMM 在凸 LP 上的实现,收敛问题

如何解决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 举报,一经查实,本站将立刻删除。