# CVXPY ECOS and SCS solver stuck in local optima for convex problem?

Consider the following optimization problem:

\begin{equation*} \begin{aligned} & \underset{z_{k, t}}{\text{maximize}} & & \sum_{k} \sum_{t} \text{concave-function}(z_{k,t}) \\ & \text{subject to} & & \text{positive lb} \leq z_{k, t} \leq \text{positive ub}, \quad \forall k, \forall t \\ &&& \sum_{k} \sum_{t} z_{k, t} \leq \text{total} \\ &&& \frac{\sum_{k}\sum_{t} z_{k, t}}{\sum_{k}\sum_{t}\text{concave-function}(z_{k, t})} \leq 10 \\ &&& z_{k,t} * 1.5 \geq z_{k, t+1} \\ &&& z_{k,t} * 0.5 \leq z_{k, t+1} \end{aligned} \end{equation*}

which should be solvable with a convex solver such as ECOS. However, when trying to solve it I end up with suboptimal solutions. Check the code below:

import numpy as np

import cvxpy as cp

bounds = [(6.274329865160556e-06, 14.919501068721157), (6.274329865160556e-06, 14.919501068721157), (6.274329865160556e-06, 14.919501068721157),
(6.274329865160556e-06, 14.919501068721157), (6.274329865160556e-06, 14.919501068721157), (6.274329865160556e-06, 14.919501068721157),
(6.274329865160556e-06, 14.919501068721157), (6.274329865160556e-06, 14.919501068721157), (6.274329865160556e-06, 14.919501068721157),
(6.274329865160556e-06, 14.919501068721157), (6.274329865160556e-06, 14.919501068721157), (6.274329865160556e-06, 14.919501068721157),
(6.274329865160556e-06, 14.919501068721157), (6.274329865160556e-06, 14.919501068721157), (6.274329865160556e-06, 14.919501068721157),
(6.274329865160556e-06, 14.919501068721157), (6.274329865160556e-06, 14.919501068721157), (6.274329865160556e-06, 14.919501068721157),
(6.274329865160556e-06, 14.919501068721157), (6.274329865160556e-06, 14.919501068721157), (6.274329865160556e-06, 14.919501068721157),
(6.274329865160556e-06, 14.919501068721157), (6.274329865160556e-06, 14.919501068721157), (6.274329865160556e-06, 14.919501068721157),
(6.274329865160556e-06, 14.919501068721157), (0.0005248360152365689, 27.172072902944304), (0.0005248360152365689, 27.172072902944304),
(0.0005248360152365689, 27.172072902944304), (0.0005248360152365689, 27.172072902944304), (0.0005248360152365689, 27.172072902944304),
(0.0005248360152365689, 27.172072902944304), (0.0005248360152365689, 27.172072902944304), (0.0005248360152365689, 27.172072902944304),
(0.0005248360152365689, 27.172072902944304), (0.0005248360152365689, 27.172072902944304), (0.0005248360152365689, 27.172072902944304),
(0.0005248360152365689, 27.172072902944304), (0.0005248360152365689, 27.172072902944304), (0.0005248360152365689, 27.172072902944304),
(0.0005248360152365689, 27.172072902944304), (0.0005248360152365689, 27.172072902944304), (0.0005248360152365689, 27.172072902944304),
(0.0005248360152365689, 27.172072902944304), (0.0005248360152365689, 27.172072902944304), (0.0005248360152365689, 27.172072902944304),
(0.0005248360152365689, 27.172072902944304), (0.0005248360152365689, 27.172072902944304), (0.0005248360152365689, 27.172072902944304),
(0.0005248360152365689, 27.172072902944304), (0.0005248360152365689, 27.172072902944304), (0.00014238207960809104, 60.130237219435635),
(0.00014238207960809104, 60.130237219435635), (0.00014238207960809104, 60.130237219435635), (0.00014238207960809104, 60.130237219435635),
(0.00014238207960809104, 60.130237219435635), (0.00014238207960809104, 60.130237219435635), (0.00014238207960809104, 60.130237219435635),
(0.00014238207960809104, 60.130237219435635), (0.00014238207960809104, 60.130237219435635), (0.00014238207960809104, 60.130237219435635),
(0.00014238207960809104, 60.130237219435635), (0.00014238207960809104, 60.130237219435635), (0.00014238207960809104, 60.130237219435635),
(0.00014238207960809104, 60.130237219435635), (0.00014238207960809104, 60.130237219435635), (0.00014238207960809104, 60.130237219435635),
(0.00014238207960809104, 60.130237219435635), (0.00014238207960809104, 60.130237219435635), (0.00014238207960809104, 60.130237219435635),
(0.00014238207960809104, 60.130237219435635), (0.00014238207960809104, 60.130237219435635), (0.00014238207960809104, 60.130237219435635),
(0.00014238207960809104, 60.130237219435635), (0.00014238207960809104, 60.130237219435635), (0.00014238207960809104, 60.130237219435635)]

maxcovers = {'x1': 2.5, 'x2': 2.5, 'x3': 2.5}
maxjumps = {'x1': 2.5, 'x2': 2.5, 'x3': 2.5}

paramvaluemappings = {"exponents": {0: np.full((1, 31), 0.034835),
1: np.full((1, 31), 0.039884),
2: np.full((1, 31), 0.083774)},
"timevariables": {0: np.full((1, 31), -1.206971),
1: np.full((1, 31), 0.36720),
2: np.full((1, 31), 1.0164)},
"lag": {0: np.zeros((1)),
1: np.zeros((1)),
2: np.zeros((1))},
"intercept": np.full((1, 31), 0.1899)}

totalvarsum = 2194

x0 = [0.01000627432986516, 0.01000627432986516, 0.01000627432986516, 0.01000627432986516,
0.01000627432986516, 0.01000627432986516, 0.01000627432986516, 0.01000627432986516,
0.01000627432986516, 0.01000627432986516, 0.01000627432986516, 0.01000627432986516,
0.01000627432986516, 0.01000627432986516, 0.01000627432986516, 0.01000627432986516,
0.01000627432986516, 0.01000627432986516, 0.01000627432986516, 0.01000627432986516,
0.01000627432986516, 0.01000627432986516, 0.01000627432986516, 0.01000627432986516,
0.01000627432986516, 0.01052483601523657, 0.01052483601523657, 0.01052483601523657,
0.01052483601523657, 0.01052483601523657, 0.01052483601523657, 0.01052483601523657,
0.01052483601523657, 0.01052483601523657, 0.01052483601523657, 0.01052483601523657,
0.01052483601523657, 0.01052483601523657, 0.01052483601523657, 0.01052483601523657,
0.01052483601523657, 0.01052483601523657, 0.01052483601523657, 0.01052483601523657,
0.01052483601523657, 0.01052483601523657, 0.01052483601523657, 0.01052483601523657,
0.01052483601523657, 0.01052483601523657, 0.010142382079608091, 0.010142382079608091,
0.010142382079608091, 0.010142382079608091, 0.010142382079608091, 0.010142382079608091,
0.010142382079608091, 0.010142382079608091, 0.010142382079608091, 0.010142382079608091,
0.010142382079608091, 0.010142382079608091, 0.010142382079608091, 0.010142382079608091,
0.010142382079608091, 0.010142382079608091, 0.010142382079608091, 0.010142382079608091,
0.010142382079608091, 0.010142382079608091, 0.010142382079608091, 0.010142382079608091,
0.010142382079608091, 0.010142382079608091, 0.010142382079608091]

maxjumps = {"x1": 1.5,
"x2": 1.5,
"x3": 1.5}
maxcovers = {"x1": 1.5,
"x2": 1.5,
"x3": 2.5}

variables = cp.Variable(len(x0))
constraints = []
for lowerbound, upperbound, sumvariables in zip([x[0] for x in bounds],
[x[1] for x in bounds], variables):
constraints += [lowerbound <= sumvariables, sumvariables <= upperbound]

total_varsum_constraint = [cp.sum(variables) <= totalvarsum]
constraints += total_varsum_constraint

def newobjective2(x, *args):
paramvaluemappingsordered = args[0]
horizon = args[3]
channels = args[4]

lookbacklag = 4

previousxs = {0: np.array([3.06926, 6.40988, 12.00676, 11.4556, 32.0006, 12.0006]),
1: np.array([3.06926, 6.40988, 12.00676, 11.4556, 32.0006, 12.0006]),
2: np.array([3.06926, 6.40988, 12.00676, 11.4556, 32.0006, 12.0006])}
summed = 0
prod = 0
for idx, channel in enumerate(channels):
exponententries = paramvaluemappingsordered['exponents'][idx]
expenditurechannel = np.array(list([*previousxs[idx], *x[idx * horizon: idx * horizon + horizon]]))
channelentry = np.power(expenditurechannel, exponententries)
prod += channelentry

intercept = paramvaluemappingsordered[
'intercept'
]

prod += intercept
summed = prod

summed = summed[:, lookbacklag:]
meanvalue = np.mean(summed, axis=0)

return -np.sum(meanvalue)

breakpoints = {channel: idx * 25 for idx, channel in enumerate(["x1", "x2"])}

horizon = 25
channels = ["x1", "x2", "x3"]
minspendtotalconstraints = {"x1": 0,
"x2": 0,
"x3": 0}

constraints = []
for lowerbound, upperbound, sumvariables in zip([x[0] for x in bounds],
[x[1] for x in bounds], variables):
constraints += [lowerbound <= sumvariables, sumvariables <= upperbound]

total_varsum_constraint = [cp.sum(variables) <= totalvarsum]
constraints += total_varsum_constraint

breakpoints = {channel: idx * horizon for idx, channel in enumerate(channels)}

for channel, cutoff in breakpoints.items():
fromx = cutoff
tox = cutoff + horizon

maxjump = maxjumps[channel]
maxcover = maxcovers[channel]

minspendchannel = minspendtotalconstraints[channel]
minspendconstrainttotal = [cp.sum(variables[fromx: tox]) >= minspendchannel]
constraints += minspendconstrainttotal

for idx in range(fromx, tox - 1):
maxjumpconstraint_ = [variables[idx] * maxjump <= variables[idx + 1]]
constraints += maxjumpconstraint_

for idx in range(fromx, tox - 1):
maxcoverconstraint_ = [variables[idx + 1] * maxcover >= variables[idx]]
constraints += maxcoverconstraint_

targetvalue = 10.1

def targetpenalty(x, target, args_obj):
resp = -newobjective2(x, *args_obj)
return (resp * target) - cp.sum(x)

constraints += [targetpenalty(variables, targetvalue, (paramvaluemappings, None, None, 25, ["x1", "x2", "x3"])) >= 0]

obj = cp.Minimize(newobjective2(variables, *(paramvaluemappings, None, None, 25, ["x1", "x2", "x3"])))

prob = cp.Problem(obj, constraints)

solved = prob.solve(solver="ECOS", verbose=True)

solution = prob.solution.primal_vars[1]

sumvariables = np.sum(solution)

result = -newobjective2(solution, *(paramvaluemappings, None, None, 25, ["x1", "x2", "x3"]))

print(sumvariables / result)
print(solved)


Since it has not reached the limit of the 3rd constraint it should increase the values for x, it is easy to notice that the solutions is not optimal. Now why is that? it is a convex problem..

• Convex solvers are never stuck in local optima. Feb 7 at 16:53
• What is the "concave-function"? Is it always nonnegative (if not, the problem may not be convex)? Did you move it to the Right Hand Side of the constraint? Perhaps you can show us the solver and CVXPY output? Does CVXPY claim the problem was solved to optimality? DId you use DCP rules, as opposed to DCCp or some other "non-convex" rules which don't guarantee global optimality? Feb 7 at 16:58
• it is indeed nonnegative, it is moved to the right hand in the constraint in the code, yes. if i get rid of the constraint governing how much z_kt can increase from a timestep to another it arrives at something that seems reasonable, if not, it keeps increasing the value from some startvalue i dont understand where it gets it from, since i have not constrained z_k,t at t=0 except for the bounds it should start of with the "optimal" value as when i get rid of this constraint. If i increase it to something > 2 it does not find a solution saying infeasible, opt_val =inf which i dont understand Feb 7 at 17:46
• Mosek is likely to be more robust than the other optimizers. Feb 8 at 5:32
• I'm not sure what you mean by "3rd constraint" and what exactly the issue is, but after running your code with two different solvers I always get the same very good solution with minimal violations of order 1e-9, ie. everything looks numerically as good as one could possibly expect. Maybe you are modeling something else than you thought. Feb 8 at 8:16