ADMM implementation on a convex LP, convergence issue.

[Hassan Yazdani]

Hi. I hope you all are in good health. 

I have coded a simple LP on cvxpy, and it has an optimal answer. By relaxing one constraint, (p_p[0] + p_p[1] == 0), it can be solved as two separate problems, and I have done it using ADMM. this is done in a for loop, solving the two problems (prob_0 and prob_1) in parallel and updating the dual variable (mu_bar) in each iteration. For some reason that I haven't figured out yet, the ADMM implementation keeps oscillating on two values and I don't get a convergence. You can see the code in the colab link below:

https://colab.research.google.com/drive/1WICekkYukrrezAVb1Yz-IP2hBfX3Qnpd?usp=sharing

Any insight on how to tackle the issue is highly appreciated. Thanks. 

The ADMM code:

%%%%%%%%%%%%%%%%%%%%%%% 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.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, p), 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, :]]).T ,

                    ]

    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])

    

# plor 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()