Cvxpy: Randomness in constraints/variables order (bad for parametric problems)

[bstellato]

The order of variables and constraints passed to the solver is not always preserved.

Every time a problem is built constraints/variables can appear in different orders which correspond to data matrices with permuted rows/columns. This makes it harder to detect if OSQP (or any other solver) problem updates can be applied correctly.

_One idea could be to force the variables and constraints IDs to be in increasing/decreasing order so that the problem rows/columns are not permuted._ What do you think?

As @gbanjac pointed out, removing randomness from the data building process would make warm starting much more efficient. At the moment, OSQP needs to refactor the KKT matrix whenever the ordering of variables or constraints is not preserved.

Related to #620.

[SteveDiamond]

This should definitely be fixed. It came up in the old cvxpy, and was fixed there. I took exactly the approach you propose, where I sorted the variables and constraints. Probably we're getting the issue because cvxpy 1.0 creates a lot of new variables as it rewrites problems. Another solution would be to have some sort of hash based on the actual content of the expressions, which could be used for consistent ordering.

[dave31415]

Is this the cause of issue #614 Non-deterministic output. I've asked for that issue to be reopened as it is definitely a CVXPY issue not just an ECOS issue. See updates.

[dave31415]

It would be useful I think to guarantee that the output is deterministic. Perhaps we could add some test cases to ensure that. Otherwise it's hard to use this in production because, in the enterprise, it's often required that you can reproduce any issue exactly. And it makes testing everything else far easier.

[nepslor]

Experienced the same problem. It took me ages to find out it was a non deterministic solver behavior, since I would expect no source of randomness in a solver.

[bstellato]

Can this be relevant for Python 3.6+? The insertion order in dictionaries seems to be preserved https://stackoverflow.com/questions/39980323/are-dictionaries-ordered-in-python-3-6

[SteveDiamond]

It does seem to happen in Python 3.6 somehow. We've had reports from the class. #664 should remove any non-determinism in solving the same problem repeatedly within a piece of code, but there is still non-determinism between code runs.

[bstellato]

Could it be fixed simply by adjusting the order of what is returned by problem.variables() and problem.constraints? Maybe I am missing something.

[bstellato]

Are there any updates on this? Several people are having the same problem.

I would be happy to help but I am not sure where exactly the nondeterminism arises in the chain. I suspect (definitely not sure) in the higher level reductions in the chain and not in the matrix stuffing. @SteveDiamond do you think it is something easily fixable?

[bstellato]

Here is a mwe where I define a mixed-integer control problem and the behavior is non deterministic.
This happens if I just extract the data twice.

This makes the exact same code return different optimal solutions at each call! (If I then solve it with GUROBI)

Here is the code. Below there is the output.

EDIT: No need to solve the problem again. Just running get_problem_data twice creates the behavior within the same script. I believe it is a big reliability issue.

Code

import cvxpy as cp
import numpy as np
import scipy.sparse as spa


'''
Define problem
'''
T = 20
tau = 1.0
alpha = 6.7 * 1e-04
beta = 0.2
gamma = 0.08
E_min = 5.2
E_max = 10.2
P_max = 1.2
n_switch = 5

# Initial values
E_init = cp.Parameter(name="E_init")
z_init = cp.Parameter(name="z_init")
s_init = cp.Parameter(name="s_init")

# Load profile
P_load = cp.Parameter(T, name="P_load")

# Past d values
past_d = cp.Parameter(T, "past_d")

# Variables
E = cp.Variable(T, name="E")  # Capacitor
P = cp.Variable(T, name="P")  # Cell power
s = cp.Variable(T + 1, name="s")
z = cp.Variable(T, name="z")
w = cp.Variable(T, name="w")
d = cp.Variable(T, name="d", boolean=True)

# Constraints
bounds = []
bounds += [E_min <= E, E <= E_max]
bounds += [0 <= P, P <= z * P_max]
bounds += [-1 <= w, w <= 1]
bounds += [0.0 <= z, z <= 1]

# Capacitor dynamics
capacitor_dynamics = []
for t in range(T - 1):
    capacitor_dynamics += [E[t + 1] == E[t] + tau * (P[t] - P_load[t])]

# Switch positions
switch_positions = [z[t + 1] == z[t] + w[t] for t in range(T - 1)]

switch_accumulation = []
for t in range(T):
    switch_accumulation += [s[t + 1] == s[t] + d[t] - past_d[t]]

# Number of switchings
number_switchings = [s <= n_switch]

# Logical relationships
G = np.array([[1.0, 0.0, -1.0],
              [-1.0, 0.0, -1.0],
              [1.0, 2.0, 2.0],
              [-1.0, -2.0, 2.0]])
g = np.array([0.0, 0.0, 3.0, 1.0])

logics = []
for t in range(T):
    logics += [G[:, 0] * w[t] + G[:, 1] * z[t] + G[:, 2] * d[t] <= g]

# Initial values (parameters)
initial_values = [E[0] == E_init, z[0] == z_init, s[0] == s_init]

constraints = (bounds + capacitor_dynamics + initial_values
               + switch_positions + switch_accumulation
               + number_switchings + logics)

# Objective
cost = 0
for t in range(T - 1):
    cost += alpha * (P[t]) ** 2 + beta * P[t] + gamma * z[t]

objective = cp.Minimize(cost)

problem = cp.Problem(objective, constraints)

'''
Define parameters
'''
params = {}
params["E_init"] = 6.11309086944192
params["z_init"] = 0.0
params["s_init"] = 0.0
params["past_d"] = np.zeros(T)
params["P_load"] = np.array([0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
                             0.0, 0.01161157, 0.02576362, 0.03991567,
                             0.05476047, 0.07064437, 0.08652828,
                             0.10241218, 0.12006155, 0.14035911,
                             0.16065667, 0.18095423, 0.18386594, 0.16069887])

'''
Nondeterministic behavior
'''
solver = cp.GUROBI

# Populate parameters and get_data
for p in problem.parameters():
    p.value = params[p.name()]

data_first = problem.get_problem_data(solver)[0]
data_second = problem.get_problem_data(solver)[0]

print("Difference returned data\n---")
for k in data_first.keys():
    if spa.isspmatrix(data_first[k]):
        difference = (data_first[k] - data_second[k]).todense()
    elif not isinstance(data_first[k], dict):
        difference = np.array(data_first[k]) - np.array(data_second[k])
    diff = np.linalg.norm(difference)

    print("Difference in %s: %.2e" % (k, diff))

Output

Difference returned data
---
Difference in P: 5.69e-03
Difference in q: 2.84e-03
Difference in A: 0.00e+00
Difference in b: 0.00e+00
Difference in F: 0.00e+00
Difference in G: 0.00e+00
Difference in bool_vars_idx: 0.00e+00
Difference in int_vars_idx: 0.00e+00
Difference in n_var: 0.00e+00
Difference in n_eq: 0.00e+00
Difference in n_ineq: 0.00e+00

[bstellato]

This might be of interest: https://stackoverflow.com/a/54545521/9753175

Dictionaries preserve insertion order from Python 3.6 onwards but sets do not.

[SteveDiamond]

I investigated the example above, and the problem is not with variable orderings but with the actual data changing. Namely, the first call to get_problem_data returns P and q with different entries, not just rearranged entries. I don't know why this would happen.

[SteveDiamond]

I fixed the example above on master. I'll close this after merging #698.

[bstellato]

Thank you for looking into this! Just a curiosity, does this fix affect/solve https://github.com/cvxgrp/cvxpy/issues/614?

[SteveDiamond]

Yes, the code in #614 works for me.

[hroskes]

Hi, I just had a similar problem and came across this issue in my search. I upgraded to master and that fixed the problem. Any chance of a new release relatively soon so that other people don't waste more time on a problem that's already solved in master?

[SteveDiamond]

I pushed a new release. I hadn't realized it had been so long. It will be out on PyPi and the cvxgrp conda channel pretty soon.

[hroskes]

Thanks!