How to do 'if x, then y == z' as a constraint?

[Matthew Davis]

Background:

This page explains how to implement

if x:
   y > z

As a constraint, where x,y and z are cvxpy variables. The solution is to add the constraint

(1-x) * large_num > z - y

Where large_num is far larger than z and y.
The way it works is that

• if x is false (0), then we get "large_num > z-y", which is always true, so the constraint is redundant.
• if x is true (1), then we get 0 > z-y   which is equivalent to y > z

It's a bit confusing, but it is valid.

My problem:

I want to do exactly that, but instead of an inequality (y > z), I want to do an equality (y == z)

So I want to write a constraint to implement:

if x:
   y == z

How can I do that?

I am trying to implement:

if x:
    y >= z
    AND
    y <= z

Since (y <=z) and (y >=z) are both true only if y == z, which is what I want.

To write this properly, I have:

x = cvxpy.Int(name='x')
y = cvxpy.Int(name='y')
z = cvxpy.Int(name='z')

large_num = 99999 # far larger than any realistic value of y or z

constraints = []

constraints.append((1-x) * large_num >= y - z)
constraints.append((1-x) * -1 * large_num <= y - z)

...

When I try this, the solver doesn't throw an error, and claims to have found the 'optimal' solution, but the constraint 'if x then y == z' is violated by the solution.

More specifically, it works when this is the only constraint, but when I have two independent constraints in this way (where it happens that at most one should be active at a time) then at least one of those constraints is violated. But the solver doesn't complain.

What is wrong with my approach? Is there a simpler way?

(I have tried with x as a bool too. Same result.)

[Steven Diamond]

Your code/approach looks correct. Are you sure the status isn't "OPTIMAL_INACCURATE"? Are you using the ECOS_BB solver (set verbose=True)? Try using a different MIP solver. To say any more I need a self contained snippet that replicates the problem.

[Matthew Davis]

Here's a minimal working example.

I've included print statements to print out the states which violate the constraint, and the constraint which is violated.

Interestingly, problem.status is 'optimal'. But 'problem.status == cvxpy.OPTIMAL' returns false. Why is that?
The verbose output says 'PRIMAL INFEASIBLE'.
So is that optimal or not?

import numpy
import cvxpy
from cvxopt.modeling import op
from cvxopt.modeling import variable, op, max, sum

num_states = 10 # each cvxpy variable is an array this long

min_on_time = 3 # if states turns from 0 to 1, it must remain 1 for at least 3 states

states = cvxpy.Bool(num_states,name='states')

time_remaining = cvxpy.Variable(num_states,name='time_remaining')

constr = [] # constraints list

# just forcing states to do something interesting,
# to get the behaviour I'm talking about
#for (i,s) in enumerate(states):
#    if (i%5) <= 2:
#        constr.append(s == 1)
#    else:
#        constr.append(s == 0)

# setup the first period
constr.append(time_remaining[0] == min_on_time)

# now is the tricky bit.
# I want to implement:
# if states[i-1]:
#    time_remaining[i] = time_remaining[i-1] - 1
# else:
#    time_remaining[i] = min_on_time * states[i]

large_num = num_states * 1e3

for i in range(num_states):
    # if states[i-1] then time_remaining[i] = time_remaining[i-1] - solve_period_len
    constr.append((1-states[i-1]) * large_num   >= time_remaining[i-1] - 1   - time_remaining[i])
    constr.append((1-states[i-1])*-1*large_num  <= time_remaining[i-1] - 1  - time_remaining[i])

    # if not states[i-1] then time_remaining[i] = min_on_time * states[i]
    constr.append(states[i-1] * large_num    >= min_on_time * states[i] - time_remaining[i])
    constr.append(states[i-1] * -1*large_num <= min_on_time * states[i] - time_remaining[i])

for i in range(num_states):
     constr.append(time_remaining[i] <= states[i] * large_num)


# setup some objective, I don't care what
obj = []

for (i,s) in enumerate(states):
    if i%2:
        obj.append(cvxpy.Maximize(s))
    else:
        obj.append(cvxpy.Maximize(-1*s))


problem = cvxpy.Problem(sum(obj), constr)

ret = problem.solve(verbose=True)

print('problem.status == %s' % problem.status)
print('problem.status == cvxpy.OPTIMAL :  %s' % problem.status == cvxpy.OPTIMAL)


print('states:         [%s]' % ','.join(['%2d' % round(x.value) for x in states]))
print('time_remaining: [%s]' % ','.join(['%2d' % round(x.value) for x in time_remaining]))

for i in range(num_states):
    if round(states[i-1].value):
        if round(time_remaining[i].value) != round(time_remaining[i-1].value) - 1:
            print('invalid for i=%d' % i)
            print('    constraint violated: if states[%d], then time_remaining[%d]==time_remaining[%d]-1' % (i,i,i-1))
            print('    where:')
            print('        states[%d]=%f' % (i,states[i].value))
            print('        time_remaining[%d]=%f' % (i,time_remaining[i].value))
            print('        time_remaining[%d]=%f' % (i-1,time_remaining[i-1].value))
    else:
        if round(time_remaining[i].value) != (min_on_time * round(states[i].value)):
            print('invalid for i=%d' % i)
            print('    constraint violated: time_remaining[%d] = min_on_time * states[%d]' % (i,i))
            print('    where:')
            print('        time_remaining[%d]=%f' % (i,time_remaining[i].value))
            print('        min_on_time = %f' % min_on_time)
            print('        states[%d]=%f' % (i,states[i].value))
            print('        min_on_time * states[%d] = %f' % (i,min_on_time*states[i].value))

The output (only including the verbose output from the last iteration) is

ECOS 2.0.4 - (C) embotech GmbH, Zurich Switzerland, 2012-15. Web: www.embotech.com/ECOS

It     pcost       dcost      gap   pres   dres    k/t    mu     step   sigma     IR    |   BT
 0  -1.442e-02  -2.005e+03  +4e+03  1e-01  4e-01  1e+00  5e+01    ---    ---    1  1  - |  -  -
 1  +1.303e-01  -3.231e+03  +3e+03  3e-01  3e-01  2e+01  5e+01  0.2961  7e-01   1  0  0 |  0  0
 2  +2.458e-02  -5.508e+02  +7e+02  5e-02  4e-02  1e+01  9e+00  0.8410  5e-02   1  0  0 |  0  0
 3  +2.324e-03  -3.061e+01  +4e+01  3e-03  2e-03  9e-01  5e-01  0.9609  2e-02   1  0  0 |  0  0
 4  +5.093e-05  -5.525e-01  +6e-01  5e-05  4e-05  2e-02  9e-03  0.9821  1e-04   1  0  0 |  0  0
 5  +3.709e-05  -2.284e-01  +2e-01  2e-05  1e-05  8e-03  4e-03  0.6630  7e-02   1  0  0 |  0  0
 6  +2.897e-04  -1.789e+00  +1e-01  3e-04  6e-05  2e-01  2e-03  0.9104  4e-01   1  0  0 |  0  0
 7  +7.147e-05  +6.137e-01  +3e-03  6e-05  8e-06  9e-01  4e-05  0.9787  1e-04   1  1  1 |  0  0
 8  +8.071e-05  +7.948e+01  +3e-05  5e-05  2e-05  8e+01  4e-07  0.9890  1e-04   1  0  0 |  0  0
 9  +8.003e-05  +7.164e+03  +3e-07  5e-05  4e-05  7e+03  5e-09  0.9890  1e-04   1  0  0 |  0  0

PRIMAL INFEASIBLE (within feastol=3.4e-10).
Runtime: 0.000364 seconds.

problem.status == optimal
False
states:         [ 1, 1, 1, 1, 0, 1, 1, 1, 0, 1]
time_remaining: [ 3, 2, 1, 0,-1, 3, 2, 1, 0, 4]
invalid for i=9
    constraint violated: time_remaining[9] = min_on_time * states[9]
    where:
        time_remaining[9]=3.993470
        min_on_time = 3.000000
        states[9]=0.999999
        min_on_time * states[9] = 2.999998

states[8] == False, and states[9] == True, so time_remaining[9] should be min_on_time = 3, but it's not, time_remaining[9] == 4.
Therefore the following constraint is violated for i==9.

    constr.append(states[i-1] * -1*large_num <= min_on_time * states[i] - time_remaining[i])


I'm wondering whether the issue is something to do with floating point errors.
Some of my variables are cvxpy integers, but when I print them out with %f, they aren't integers. (e.g. 2.999998)

If I make large_num larger (equal to 1e9), the 'solution' violates the constraints multiple times.
If I make large_num smaller (equal to 100), the solver takes a long time, and is 'optimal inaccurate', but returns a valid solution. (Not sure if it's optimal or not.)

Here's some context for this minimal example:

states is an array of bools, which is num_states=10, I haven't constrained it in any way in this example.

states is representing an electrical generator, which is either on or off. This generator has a constraint, that if switched on, should remain on for at least min_on_time=3 periods.
To implement this, I have a variable my_vars.
When the generator is switched from off to on (states[i-1] ==0, states[i] == 1), time_remaining=min_on_time.
When the generator stays on (states[i-1]==states[i]==1), or transitions from on to off (states[i-1]==1,states[i]==0), time_remaining is decremented by one.
When the generator stays off (states[i-1] == states[i]==0) I don't really mind what happens to time_remaining, but ideally time_remaining[i] = 0.

I set the objective I did, just to get the states to change in a way that exhibits the behaviour I'm talking about.

[Matthew Davis]

Sorry, my last message was truncated.

The MWE is as follows

# end of MWE

[Steven Diamond]

This is very strange behavior indeed. The ECOS_BB solver finds an optimal solution, but maybe it's confused because of floating point issues. I strongly strongly recommend if you care about this problem seriously that you use a more professional grade MIP solver like Cbc or MOSEK. If you still have problems feel free to post again.

[Matthew Davis]

There was a bug in that MWE. The for loop was supposed to have an if statement to exclude i==0, (since the contents of the for loop index i-1).

When I fixed that bug, it claims the problem is infeasible.

I've tried rephrasing this problem in several different ways.
Each time I think I've got a problem which is valid, the solver says it's infeasible.

I'm wondering. If:

• I don't specify which solver to use
• the problem is DCP
  
• the constraints are all linear
• a solution does exist

will the solver definitely return a solution?

Will the solver ever say the problem is infeasible if a solution exists?

I have looked into CBC. I must say, the installation of CBC is incredibly confusing.

I installed coin with

apt-get install coinor-cbc

But then when I say 

ret = problem.solve(solver=cvxpy.CBC)

it raises an error, saying that CBC is not installed.

Do I need to use the library cylp?
I tried installing the python module 'cylp' with pip (in a virtual env). That install fails because the os environment variable COIN_INSTALL_DIR is not set.

What's that variable supposed to be set to if I installed COIN with apt-get? (I also posted that question on the cylp github page)

I then tried compiling cylp from source, which took a while because some dependencies weren't automatically installed.

Then once I'd finished that, I tried running my original script again, and it says CBC is not installed.

When I try `import cylp` at the top of the script, it throws an error because cylp is supposedly not installed.

I'm quite lost at the moment. I don't even know if I should be asking for help from coin, cylp or cvxpy. 
How am I supposed to solve problems in cvxpy with CBC?

[Steven Diamond]

The solver should never return infeasible if a solution exists. Surely you can remove constraints until it works or construct a problem that you 100% know has a solution?

This script shows how to install Cbc on ubuntu:

https://github.com/cvxgrp/cvxpy/blob/master/continuous_integration/install.sh

[Matthew Davis]

I tried removing constraints until it works.
What I found is that there are no constraints which are invalid.

e.g. consider the constraints x<=0 and x>=1.

Neither constraint is wrong, but together they create an infeasible solution, and neither constraint is more 'wrong' than the other.

I suspect that's what's happening, but with a far more complex space. I've tried looking for such conflicts, but to no 

The problem I'm trying to solve definitely has a solution. So I probably have just written one of the constraints slightly wrong.

Hmm. I'll see if I can write an even more minimal working example.

As far as Cbc goes, I just ran that script. (Modified a bit, because some directory paths were hard coded to the user, and because there was no env to 'deactivate')

Now I have successfully installed cvxpy, cvxopt and cylp in a virtual env, and I have imported them too.

But I still get the error "The solver CBC is not installed". Is there some trick I'm missing?

I had a setup which works with the normal solver (when some of my constraints are commented out). Then I imported cylp, and called

problem.solve(solver=cvxpy.CBC)

Is there anything else I need to do to use CBC as a solver?

[Steven Diamond]

The test for whether Cbc is installed is running this code without error: "from cylp.cy import CyClpSimplex"

If you can do that in your Python environment, cvxpy should be able to call Cbc. If it can't then cvxpy and cylp are installed in different environments.

[Matthew Davis]

Ah, I can't do that.

I'm a bit new to virtual envs, so the problem is probably there.

Here's the output of my terminal when I:

• activate the virtual-env
• install cylp, cvxpy and cvxopt with pip in the virtual-env
• run the command you suggested inside a python interpreter

matt@machine:~/muckaround/opt$ . ./worker-env/bin/activate

(worker-env) matt@ip-10-61-176-23:~/muckaround/opt$ ./worker-env/bin/pip install cylp

Requirement already satisfied: cylp in ./worker-env/lib/python2.7/site-packages

Requirement already satisfied: numpy>=1.5.0 in ./worker-env/lib/python2.7/site-packages (from cylp)

Requirement already satisfied: scipy>=0.10.0 in ./worker-env/lib/python2.7/site-packages (from cylp)

(worker-env) matt@machine:~/muckaround/opt$ ./worker-env/bin/pip install cvxopt

Requirement already satisfied: cvxopt in ./worker-env/lib/python2.7/site-packages

(worker-env) matt@machine:~/muckaround/opt$ ./worker-env/bin/pip install cvxpy

Requirement already satisfied: cvxpy in ./worker-env/lib/python2.7/site-packages

Requirement already satisfied: six in ./worker-env/lib/python2.7/site-packages (from cvxpy)

Requirement already satisfied: toolz in ./worker-env/lib/python2.7/site-packages (from cvxpy)

Requirement already satisfied: multiprocess in ./worker-env/lib/python2.7/site-packages (from cvxpy)

Requirement already satisfied: scs>=1.1.3 in ./worker-env/lib/python2.7/site-packages (from cvxpy)

Requirement already satisfied: ecos>=2 in ./worker-env/lib/python2.7/site-packages (from cvxpy)

Requirement already satisfied: numpy>=1.9 in ./worker-env/lib/python2.7/site-packages (from cvxpy)

Requirement already satisfied: fastcache in ./worker-env/lib/python2.7/site-packages (from cvxpy)

Requirement already satisfied: CVXcanon>=0.0.22 in ./worker-env/lib/python2.7/site-packages (from cvxpy)

Requirement already satisfied: scipy>=0.15 in ./worker-env/lib/python2.7/site-packages (from cvxpy)

Requirement already satisfied: dill>=0.2.6 in ./worker-env/lib/python2.7/site-packages (from multiprocess->cvxpy)

(worker-env) matt@machine:~/muckaround/opt$ python

Python 2.7.12 (default, Nov 19 2016, 06:48:10)

[GCC 5.4.0 20160609] on linux2

Type "help", "copyright", "credits" or "license" for more information.

>>> from cylp.cy import CyClpSimplex

Traceback (most recent call last):

  File "<stdin>", line 1, in <module>

  File "/home/matt/muckaround/opt/worker-env/local/lib/python2.7/site-packages/cylp/cy/__init__.py", line 1, in <module>

    from CyCoinIndexedVector import CyCoinIndexedVector

ImportError: libCoinUtils.so.0: cannot open shared object file: No such file or directory

>>>

Is that how I'm supposed to install these packages into my environment?
Am I supposed to install with pip before or after activating the virtual env?

[Matthew Davis]

Just to add to that,

When I try:

import cylp

That does not throw an error

[Matthew Davis]

Sorry for discussing two separate issues (if statement constraints, and installation of CBC) in the one thread. Should I split them up into separate threads?

In regards to the infeasible problem issue, I've found that the problem is feasible for the objective

cvxpy.Maximize(states[0] - states[1] + states[2] - states[3] + states[4] - states[5] ...)

(maximize every even numbered state, minimize every odd numbered state)

But it is infeasible for

cvxpy.Maximize(states[0] - states[1] - states[2] + states[3] - states[4] - states[5] ...)

(Maximize every 3rd state, minimize the others)

Why does changing the objective change the feasibility?

For completeness, here's the code I'm running now:

import numpy
import cvxpy
from cvxopt.modeling import op
from cvxopt.modeling import variable, op, max, sum

import cylp

# round up x to the next multiple of y
def round_up(x,y):
    if x%y == 0:
        result = x
    else:
        result = x - (x%y) + y

    #print('round(%f,%f) == %f' % (x,y,result))
    return result


def my_solve(solve_period_len, min_on_time, num_periods, prev_state, prev_on_time, obj1, obj2, tol=0,verbose=False):

    states = cvxpy.Bool(num_periods,name='states')

    time_remaining = cvxpy.Variable(num_periods,name='time_remaining')

    constr = [] # constraints list

    large_num = num_periods * 1e2

    # setup the first period

    if prev_state:
        constr.append(time_remaining[0] == min_on_time - prev_on_time)
    else:
        constr.append(time_remaining[0] == min_on_time * states[0])

    for i in range(num_periods):
        if i > 0:

            # if states[i-1] then time_remaining[i] = time_remaining[i-1] - solve_period_len

            constr.append((1-states[i-1]) * large_num   >= time_remaining[i-1] - solve_period_len   - time_remaining[i])
            constr.append((1-states[i-1])*-1*large_num  <= time_remaining[i-1] - solve_period_len  - time_remaining[i])

            # if not states[i-1] then time_remaining[i] = min_on_time * states[i]
            constr.append(states[i-1] * large_num    >= min_on_time * states[i] - time_remaining[i])
            constr.append(states[i-1] * -1*large_num <= min_on_time * states[i] - time_remaining[i])

    for i in range(num_periods):

         constr.append(time_remaining[i] <= states[i] * large_num)

    # setup the objective function
    obj = cvxpy.Maximize(0)
    for (i,s) in enumerate(states):
        if (i%obj1) < obj2:
            obj += cvxpy.Maximize(states[i])
        else:
            obj += cvxpy.Maximize(-1*states[i])

    problem = cvxpy.Problem(obj, constr)

    ret = problem.solve()




    ##########
    # Everything beyond here is just printing and verifying the answer
    #########

    print('problem.status == %s' % problem.status)

    print('objectives: %s' % ','.join(['%d' % (2*(i%obj1 < obj2)-1) for i in range(num_periods)]))

    rem_calc_valid_all  = True
    sw_off_constr_valid_all = True

    for i in range(num_periods):

        row_strs = ['period:%2d' % i]

        if states[i].value == None:
            print('states[%d]==None' % i)
            return 0
        state = round(states[i].value)
        row_strs.append('state:%1d' % state)

        t_rem_i = time_remaining[i].value
        if t_rem_i != None:
            row_strs.append('t_rem:%4d' % round(t_rem_i))
        else:
            row_strs.append('t_rem:None')

        obj_i = (2*(i%obj1 < obj2)-1)

        row_strs.append('obj:%d' % obj_i)

        if i == 0:
            prev_state_i = prev_state
        else:
            prev_state_i = states[i-1].value

        # has the t_rem decreased and increased properly?
        rem_calc_valid_i = True
        if i == 0:
            if prev_state_i:
                if round(time_remaining[i].value) != (min_on_time - prev_on_time):
                    print('HERE A')
                    rem_calc_valid_i = False
            elif round(states[0].value):
                if round(time_remaining[i].value) != min_on_time:
                    print('HERE B')
                    rem_calc_valid_i = False
            elif round(time_remaining[i].value) != 0:
                rem_calc_valid_i = False
        if i > 0:
            if round(states[i-1].value):
                if round(time_remaining[i].value) != round(time_remaining[i-1].value) - solve_period_len:
                    rem_calc_valid_i = False
            else:
                if round(states[i].value) and round(time_remaining[i].value) != (min_on_time * round(states[i].value)):
                    rem_calc_valid_i = False

        row_strs.append('rem_calc_valid: %d' % rem_calc_valid_i)

        # has the device switched off when t_rem > 0?
        sw_off_constr_valid_i = True

        if prev_state_i and (not round(states[i].value)) and (round(time_remaining[i].value) > 0):
            sw_off_constr_valid_i = False

        row_strs.append('sw_off_constr: %d' % sw_off_constr_valid_i)

        overall_valid_i = rem_calc_valid_i and sw_off_constr_valid_i

        row_strs.append('overall: %d' % overall_valid_i)

        print(' | '.join(row_strs))

        if not rem_calc_valid_i:
            rem_calc_valid_all = False
        if not sw_off_constr_valid_i:
            sw_off_constr_valid_all = False

    overall_valid = rem_calc_valid_all and sw_off_constr_valid_all

    print('\nOverall valid: %d' % overall_valid)
    return overall_valid


if __name__ == "__main__":
    #main()
    ret = my_solve(1, 3, 30, 0, 0, 3, 1)
    assert(ret)

[Steven Diamond]

You should use Anaconda, not virtualenv. Installing everything with virtualenv is a nightmare that I can't help you with. I'm happy to take another look at your problem once you've tried solving it with a robust MIP solver like Cbc, MOSEK, GUROBI, etc.