Simple constrained optimization problem with mystic gives suboptimal results and violates constraints

Question

I am working on a optimization problem and am wondering why the solution is so far from optimum.

As a part of it i made the objective function much easier so it is easy to analytically tell where the solver should allocate more resources: in this simplified version it is easy to tell that i should maximize w.r.t constraints inv2 for days near the end of the month since the coefficient for inv2 is higher than inv1 and the day coefficient is positive.

The solver does not solve this however. It also seems strange that it does not use the whole resource = 100000 thus violates this constraint(last one).

I think i have made a small but vital error in the code somewhere.

So Before digging deep into differential evolutional algorithms i thought i ask here.

Why does the solver violate the constraint and find such an suboptimal solution?

import pandas
from datetime import timedelta
from datetime import datetime
from dateutil.tz import gettz
import mystic
from mystic.symbolic import simplify
from mystic.symbolic import generate_constraint, generate_solvers
from mystic.monitors import VerboseMonitor
from mystic.symbolic import symbolic_bounds


horizon = 30
nowianatimezone = datetime.now(tz=gettz(name='Europe/Stockholm'))
datesbetweendates = pandas.date_range(start=nowianatimezone,
                                      end=nowianatimezone + timedelta(days=horizon),
                                      freq='d')
daysforcalculation = datesbetweendates.day.tolist()

def objectiveforoneday(inv1, inv2, day):
    return 1.1 * inv1 + 1.4 * inv2 + 10 * day
            
def objective(x):
    tbret = 0
    for idx, day in enumerate(daysforcalculation):
        inv1 = x[idx * 2]
        inv2 = x[(idx * 2) + 1]
        tbret += objectiveforoneday(inv1=inv1,
                                    inv2=inv2,
                                    day=day)
    return tbret

constraints = '''
x0 >= 250.0
x1 >= 100.0
x2 >= 250.0
x3 >= 100.0
x4 >= 250.0
x5 >= 100.0
x6 >= 250.0
x7 >= 100.0
x8 >= 250.0
x9 >= 100.0
x10 >= 250.0
x11 >= 100.0
x12 >= 250.0
x13 >= 100.0
x14 >= 250.0
x15 >= 100.0
x16 >= 250.0
x17 >= 100.0
x18 >= 250.0
x19 >= 100.0
x20 >= 250.0
x21 >= 100.0
x22 >= 250.0
x23 >= 100.0
x24 >= 250.0
x25 >= 100.0
x26 >= 250.0
x27 >= 100.0
x28 >= 250.0
x29 >= 100.0
x30 >= 250.0
x31 >= 100.0
x32 >= 250.0
x33 >= 100.0
x34 >= 250.0
x35 >= 100.0
x36 >= 250.0
x37 >= 100.0
x38 >= 250.0
x39 >= 100.0
x40 >= 250.0
x41 >= 100.0
x42 >= 250.0
x43 >= 100.0
x44 >= 250.0
x45 >= 100.0
x46 >= 250.0
x47 >= 100.0
x48 >= 250.0
x49 >= 100.0
x50 >= 250.0
x51 >= 100.0
x52 >= 250.0
x53 >= 100.0
x54 >= 250.0
x55 >= 100.0
x56 >= 250.0
x57 >= 100.0
x58 >= 250.0
x59 >= 100.0
x60 >= 250.0
x61 >= 100.0
x0 <= 16666.666666666668
x1 <= 16666.666666666668
x2 <= 16666.666666666668
x3 <= 16666.666666666668
x4 <= 16666.666666666668
x5 <= 16666.666666666668
x6 <= 16666.666666666668
x7 <= 16666.666666666668
x8 <= 16666.666666666668
x9 <= 16666.666666666668
x10 <= 16666.666666666668
x11 <= 16666.666666666668
x12 <= 16666.666666666668
x13 <= 16666.666666666668
x14 <= 16666.666666666668
x15 <= 16666.666666666668
x16 <= 16666.666666666668
x17 <= 16666.666666666668
x18 <= 16666.666666666668
x19 <= 16666.666666666668
x20 <= 16666.666666666668
x21 <= 16666.666666666668
x22 <= 16666.666666666668
x23 <= 16666.666666666668
x24 <= 16666.666666666668
x25 <= 16666.666666666668
x26 <= 16666.666666666668
x27 <= 16666.666666666668
x28 <= 16666.666666666668
x29 <= 16666.666666666668
x30 <= 16666.666666666668
x31 <= 16666.666666666668
x32 <= 16666.666666666668
x33 <= 16666.666666666668
x34 <= 16666.666666666668
x35 <= 16666.666666666668
x36 <= 16666.666666666668
x37 <= 16666.666666666668
x38 <= 16666.666666666668
x39 <= 16666.666666666668
x40 <= 16666.666666666668
x41 <= 16666.666666666668
x42 <= 16666.666666666668
x43 <= 16666.666666666668
x44 <= 16666.666666666668
x45 <= 16666.666666666668
x46 <= 16666.666666666668
x47 <= 16666.666666666668
x48 <= 16666.666666666668
x49 <= 16666.666666666668
x50 <= 16666.666666666668
x51 <= 16666.666666666668
x52 <= 16666.666666666668
x53 <= 16666.666666666668
x54 <= 16666.666666666668
x55 <= 16666.666666666668
x56 <= 16666.666666666668
x57 <= 16666.666666666668
x58 <= 16666.666666666668
x59 <= 16666.666666666668
x60 <= 16666.666666666668
x61 <= 16666.666666666668
x0 + x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 + x10 + x11 + x12 + x13 + x14 + x15 + x16 + x17 + x18 + x19 + x20 + x21 + x22 + x23 + x24 + x25 + x26 + x27 + x28 + x29 + x30 + x31 + x32 + x33 + x34 + x35 + x36 + x37 + x38 + x39 + x40 + x41 + x42 + x43 + x44 + x45 + x46 + x47 + x48 + x49 + x50 + x51 + x52 + x53 + x54 + x55 + x56 + x57 + x58 + x59 + x60 + x61 == 100000'''

constraintsimplified = simplify(constraints,
                                all=True)

generatedconstraints = generate_constraint(generate_solvers(constraintsimplified))

initvals = [0.0 for _ in range(len(daysforcalculation) * 2)]
stepmon=VerboseMonitor(10)
import numpy as np
cost = lambda x: -objective(x)
optimum = mystic.solvers.diffev2(cost,
                                 x0=initvals,
                                 constraints=generatedconstraints,
                                 # bounds=bounds,
                                 #full_output=True,
                                 #stepmon=stepmon,
                                 #disp=True,
                                 npop=100,
                                 maxiter=10000)

optimalvalueslst = optimum[0]

# gives us 

array([16666.66666667,   100.        ,   250.        ,   100.        ,
         250.        ,   100.        ,   250.        ,   100.        ,
         250.        ,   100.        ,   250.        ,   100.        ,
         250.        ,   100.        ,   250.        ,   100.        ,
         250.        ,   100.        ,   250.        ,   100.        ,
         250.        ,   100.        ,   250.        ,   100.        ,
         250.        ,   100.        ,   250.        ,   100.        ,
         250.        ,   100.        ,   250.        ,   100.        ,
         250.        ,   100.        ,   250.        ,   100.        ,
         250.        ,   100.        ,   250.        ,   100.        ,
         250.        ,   100.        ,   250.        ,   100.        ,
         250.        ,   100.        ,   250.        ,   100.        ,
         250.        ,   100.        ,   250.        ,   100.        ,
         250.        ,   100.        ,   250.        ,   100.        ,
         250.        ,   100.        ,   250.        ,   100.        ,
         250.        ,   100.        ])

Answer 1

I'm the mystic author. The constraints are violated because they are not orthogonal. By default, mystic solves constraints sequentially. Here's the relevant doc from generate_constraints:

Notes:
    If ``join=None``, then apply the given constraints sequentially.
    Otherwise, apply the constraints concurrently using a constraints coupler.

Warning:
    This constraint generator doesn't check for conflicts in conditions, but
    simply applies conditions in the given order. This constraint generator
    assumes that a single variable has been isolated on the left-hand side
    of each constraints equation, thus all constraints are of the form
    "x_i = f(x)". This solver picks speed over robustness, and relies on
    the user to formulate the constraints so that they do not conflict.

So, to solve the constraints in parallel you need to use join, like this:

>>> from mystic.constraints import and_
>>> generatedconstraints = generate_constraint(generate_solvers(constraintsimplified), join=and_)

Otherwise, your code should be fine. I'm running it right now, and it's pretty slow... so you might think about using a parallel map (with the map keyword). I'd also repeat the all the bounds information you have in the constraints as bounds and pass the bounds to the solver using the bounds keyword as well -- that way, the DE solver won't have to work as hard to generate solutions for the constraints that are inside the bounds.

EDIT: I came back to this, took a closer look at your code, and tried a bit harder to diagnose what's actually going on.

First, I'll take my advice of moving the equality constraint to a penalty:

>>> import pandas
>>> from datetime import timedelta
>>> from datetime import datetime
>>> from dateutil.tz import gettz
>>> import mystic
>>> from mystic.symbolic import simplify
>>> from mystic.symbolic import generate_constraint, generate_solvers
>>> from mystic.monitors import VerboseMonitor
>>> from mystic.symbolic import symbolic_bounds
>>> from mystic.constraints import and_
>>> from mystic.penalty import linear_equality, quadratic_equality
>>> 
>>> horizon = 30
>>> nowianatimezone = datetime.now(tz=gettz(name='Europe/Stockholm'))
>>> datesbetweendates = pandas.date_range(start=nowianatimezone,
...                                       end=nowianatimezone + timedelta(days=horizon),
...                                       freq='d')
>>> 
>>> daysforcalculation = datesbetweendates.day.tolist()
>>> 
>>> def objectiveforoneday(inv1, inv2, day):
...     return 1.1 * inv1 + 1.4 * inv2 + 10 * day
... 
>>> def objective(x):
...     tbret = 0
...     for idx, day in enumerate(daysforcalculation):
...         inv1 = x[idx * 2]
...         inv2 = x[(idx * 2) + 1]
...         tbret += objectiveforoneday(inv1=inv1,
...                                     inv2=inv2,
...                                     day=day)
...     return tbret
... 
>>> import numpy as np
>>> cost = lambda x: -objective(x)
>>> 
>>> bounds = [250,100]*31
>>> bounds = list(zip(bounds, [16666.666666666668]*62))
>>> constraints = '''
... x0 >= 250.0
... x1 >= 100.0
... x2 >= 250.0
... x3 >= 100.0
... x4 >= 250.0
... x5 >= 100.0
... x6 >= 250.0
... x7 >= 100.0
... x8 >= 250.0
... x9 >= 100.0
... x10 >= 250.0
... x11 >= 100.0
... x12 >= 250.0
... x13 >= 100.0
... x14 >= 250.0
... x15 >= 100.0
... x16 >= 250.0
... x17 >= 100.0
... x18 >= 250.0
... x19 >= 100.0
... x20 >= 250.0
... x21 >= 100.0
... x22 >= 250.0
... x23 >= 100.0
... x24 >= 250.0
... x25 >= 100.0
... x26 >= 250.0
... x27 >= 100.0
... x28 >= 250.0
... x29 >= 100.0
... x30 >= 250.0
... x31 >= 100.0
... x32 >= 250.0
... x33 >= 100.0
... x34 >= 250.0
... x35 >= 100.0
... x36 >= 250.0
... x37 >= 100.0
... x38 >= 250.0
... x39 >= 100.0
... x40 >= 250.0
... x41 >= 100.0
... x42 >= 250.0
... x43 >= 100.0
... x44 >= 250.0
... x45 >= 100.0
... x46 >= 250.0
... x47 >= 100.0
... x48 >= 250.0
... x49 >= 100.0
... x50 >= 250.0
... x51 >= 100.0
... x52 >= 250.0
... x53 >= 100.0
... x54 >= 250.0
... x55 >= 100.0
... x56 >= 250.0
... x57 >= 100.0
... x58 >= 250.0
... x59 >= 100.0
... x60 >= 250.0
... x61 >= 100.0
... x0 <= 16666.666666666668
... x1 <= 16666.666666666668
... x2 <= 16666.666666666668
... x3 <= 16666.666666666668
... x4 <= 16666.666666666668
... x5 <= 16666.666666666668
... x6 <= 16666.666666666668
... x7 <= 16666.666666666668
... x8 <= 16666.666666666668
... x9 <= 16666.666666666668
... x10 <= 16666.666666666668
... x11 <= 16666.666666666668
... x12 <= 16666.666666666668
... x13 <= 16666.666666666668
... x14 <= 16666.666666666668
... x15 <= 16666.666666666668
... x16 <= 16666.666666666668
... x17 <= 16666.666666666668
... x18 <= 16666.666666666668
... x19 <= 16666.666666666668
... x20 <= 16666.666666666668
... x21 <= 16666.666666666668
... x22 <= 16666.666666666668
... x23 <= 16666.666666666668
... x24 <= 16666.666666666668
... x25 <= 16666.666666666668
... x26 <= 16666.666666666668
... x27 <= 16666.666666666668
... x28 <= 16666.666666666668
... x29 <= 16666.666666666668
... x30 <= 16666.666666666668
... x31 <= 16666.666666666668
... x32 <= 16666.666666666668
... x33 <= 16666.666666666668
... x34 <= 16666.666666666668
... x35 <= 16666.666666666668
... x36 <= 16666.666666666668
... x37 <= 16666.666666666668
... x38 <= 16666.666666666668
... x39 <= 16666.666666666668
... x40 <= 16666.666666666668
... x41 <= 16666.666666666668
... x42 <= 16666.666666666668
... x43 <= 16666.666666666668
... x44 <= 16666.666666666668
... x45 <= 16666.666666666668
... x46 <= 16666.666666666668
... x47 <= 16666.666666666668
... x48 <= 16666.666666666668
... x49 <= 16666.666666666668
... x50 <= 16666.666666666668
... x51 <= 16666.666666666668
... x52 <= 16666.666666666668
... x53 <= 16666.666666666668
... x54 <= 16666.666666666668
... x55 <= 16666.666666666668
... x56 <= 16666.666666666668
... x57 <= 16666.666666666668
... x58 <= 16666.666666666668
... x59 <= 16666.666666666668
... x60 <= 16666.666666666668
... x61 <= 16666.666666666668
... '''
>>> constraintsimplified = simplify(constraints,
...                                 all=True)
>>> generatedconstraints = generate_constraint(generate_solvers(constraintsimplified), join=and_)
>>> penalty = linear_equality(lambda x: sum(x) - 100000)(lambda x: 0.0)
>>> 
>>> stepmon=VerboseMonitor(1)
>>> optimum = mystic.solvers.diffev2(cost, x0=bounds, constraints=generatedconstraints, penalty=penalty, bounds=bounds, itermon=stepmon, disp=True, npop=400, maxiter=10000, gtol=200)
Generation 0 has ChiSquare: 31264592.214017
Generation 1 has ChiSquare: 29356835.910028
Generation 2 has ChiSquare: 29356835.910028
...
Generation 709 has ChiSquare: -138045.279890
Generation 710 has ChiSquare: -138045.279890
Generation 711 has ChiSquare: -138045.279890
STOP("ChangeOverGeneration with {'tolerance': 0.005, 'generations': 200}")
Optimization terminated successfully.
         Current function value: -138045.279890
         Iterations: 711
         Function evaluations: 284800
>>> sum(optimum)
100003.69609971045
>>> penalty(optimum)
369.6099710447015
>>> optimum
array([  250.        ,   100.        ,  2517.65059232,   734.92889335,
         893.19873368,  1081.78661797,   250.        ,   100.        ,
         250.        ,   100.        ,   250.        ,   100.        ,
         328.77801842,   100.        ,  1674.29503076,  1310.62426433,
        1395.39744718, 13475.53118087,  1347.41553459,   524.45422785,
         250.        ,   100.        ,  2254.99005846,   100.        ,
         250.        ,   100.        ,   250.        ,  3881.16368612,
         250.        ,   100.        ,   263.1640174 ,   541.75906308,
        1526.5462552 , 13517.25758496,   271.88092606,   100.        ,
         250.        ,  2250.99663397,   670.18534755,   100.        ,
        1898.48965304,  2333.36140119,   250.        ,  4117.40368284,
         250.        ,   102.86560974,   250.        ,   100.        ,
         572.992929  ,   100.        ,   250.        ,  2586.40181947,
         447.28471254, 15940.82553379,   250.        ,   712.20175351,
         250.        ,   100.        ,   250.        , 11573.79909848,
        1772.01300409,  2084.05278786])

You might notice that I've used itermon so the Monitor is used, and I've set x0=bounds so the DE has better variety in the population when it starts. Also, I've increased the DE capacity to npop=400 and gtol=200, so the DE has a better chance of finding a solution. It does a decent job, but you'll notice that the penalty is still slightly violated. The equality constraint is a really hard one to satisfy for 62 "randomish" numbers. You can see this by just trying the sum in the constraint, and it still fails for DE.

>>> sums = '''x0 + x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 + x10 + x11 + x12 + x13 + x14 + x15 + x16 + x17 + x18 + x19 + x20 + x21 + x22 + x23 + x24 + x25 + x26 + x27 + x28 + x29 + x30 + x31 + x32 + x33 + x34 + x35 + x36 + x37 + x38 + x39 + x40 + x41 + x42 + x43 + x44 + x45 + x46 + x47 + x48 + x49 + x50 + x51 + x52 + x53 + x54 + x55 + x56 + x57 + x58 + x59 + x60 + x61 == 100000'''
>>> simplified = simplify(sums, all=True)
>>> generated = generate_constraint(generate_solvers(simplified), join=and_)
>>> 
>>> stepmon=VerboseMonitor(1)
>>> optimum = mystic.solvers.diffev2(cost, x0=bounds, constraints=generated, penalty=None, bounds=bounds, itermon=stepmon, disp=True, npop=400, maxiter=10000, gtol=200)
Generation 0 has ChiSquare: inf
Generation 1 has ChiSquare: inf
Generation 2 has ChiSquare: inf
...
Generation 198 has ChiSquare: inf
Generation 199 has ChiSquare: inf
Generation 200 has ChiSquare: inf
STOP("ChangeOverGeneration with {'tolerance': 0.005, 'generations': 200}")
Optimization terminated successfully.
         Current function value: inf
         Iterations: 200
         Function evaluations: 80400

Basically, it fails to find a solution to the one constraint. You will get an inf (i.e. constraints are violated) even when the sum is 99999.999999999998. Like I said, the equality constraint here is a difficult one for DE to generate... so it needs help. Help can come in terms of a penalty, or potentially, by only working over integers (as that seems like it could provide a decent solution). Or, we might try a different solver than DE.

By the way, you don't have to specify the bounds in the constraints, there are other functions in mystic that can do that. Your equality sum constraint is also available in mystic. It can be efficient to build them all in one go like you did, but maybe not as convenient. Here's an alternate version of the sum constraint and the bounds constraint. I'm also going to bring in some integer constraints to help us.

>>> from mystic.constraints import normalized, boundsconstrain
>>> normed = normalized(10000)(lambda x:x)
>>> bound = boundsconstrain(*zip(*bounds))

In the above, normed is the same as the constraint built from the equality in sums, and boundsconstrain is the same as the constraint built from all of the bounds inequalities in constraints. I won't use these, but I wanted to point them out.

So, now I'll try constraining to integers, and penalizing against the exactly matching the sum. I'm also using the tightrange keyword (it does the same thing as the boundsconstrain constraint).

>>> round = lambda x: np.round(x, 0)
>>> stepmon=VerboseMonitor(1)
>>> optimum = mystic.solvers.diffev2(cost, x0=bounds, constraints=round, penalty=penalty, bounds=bounds, itermon=stepmon, disp=True, npop=400, maxiter=10000, gtol=200, tightrange=True)
Generation 0 has ChiSquare: 31985479.600000
Generation 1 has ChiSquare: 4823431.300000
Generation 2 has ChiSquare: -125692.500000
...
Generation 384 has ChiSquare: -138311.800000
Generation 385 has ChiSquare: -138311.800000
Generation 386 has ChiSquare: -138311.800000
STOP("ChangeOverGeneration with {'tolerance': 0.005, 'generations': 200}")
Optimization terminated successfully.
         Current function value: -138311.800000
         Iterations: 386
         Function evaluations: 154800
>>> optimum
array([  354.,  2355.,   250.,   149.,   250.,   100.,   265.,  9403.,
         250.,  3730.,   250.,  2700.,  1214.,   100.,   250.,  1158.,
         542.,   757.,   257.,  3409.,  1640.,   100.,  1653.,  7921.,
        3453.,   100.,   553.,  1607.,   250.,   100.,  1865.,  8483.,
         263.,   584.,   250.,  2632.,  2680.,   450.,   250.,   341.,
         622., 11941.,   250.,   191.,   250.,  1245.,  1192.,  3421.,
         250.,   305.,   265.,   100.,   250.,   159.,   468.,   990.,
         250.,  6450.,  1046.,  7088.,   250.,   100.])
>>> sum(_)
100001.0
>>> penalty(optimum)
100.0

This works pretty well, and the sum is just barely violated. So, I think with a little more tuning of the DE parameters you should be able to get it perfect. To confirm that the constraints work, we can try a different solver with the same basic constraint setup.

>>> stepmon=VerboseMonitor(1)
>>> optimum = mystic.solvers.fmin(cost, x0=[250.0]*62, constraints=round, penalty=penalty, bounds=bounds, itermon=stepmon, disp=True, maxiter=10000, ftol=1e-8,
tightrange=True)
Generation 0 has ChiSquare: 8425665.000000
Generation 1 has ChiSquare: 8341705.800000
Generation 2 has ChiSquare: 8174565.000000
...
Generation 1602 has ChiSquare: -129877.200000
Generation 1603 has ChiSquare: -129877.200000
Generation 1604 has ChiSquare: -129877.200000
Optimization terminated successfully.
         Current function value: -129877.200000
         Iterations: 1604
         Function evaluations: 3309
STOP("CandidateRelativeTolerance with {'xtol': 0.0001, 'ftol': 1e-08}")
>>> sum(optimum)
100000.0
>>> penalty(optimum)
0.0
>>> optimum
array([1427., 1410., 1144., 1479., 3051., 1607., 1248., 1002.,  840.,
       1746.,  720., 1828., 1724.,  854., 4942., 1615., 1463.,  683.,
       5441., 1590., 1226., 1489., 1637., 1595.,  773., 1646., 1366.,
       2592., 1814., 1785., 1699., 1445.,  719., 1851.,  930., 1445.,
       1034., 2065., 1728., 1502.,  729., 1900., 1644., 1905., 1759.,
       1676., 2607., 1640., 1235., 1604.,  601., 1845., 1420., 1271.,
       1628., 1736., 1794., 1614., 1471., 1542.,  462., 1762.])
>>> 

So, in this case, the constraints are perfectly satisfied, but it might be in a local minimum as I used a non-global solver. Either way, the last two runs tell me that with a bit more work tuning the solver, or manipulating how the constraints/penalties/bounds are applied, this is a hard, but solvable problem, and with tuning of the DE parameters, or potentially using one of the ensemble solvers, you can get the global optimum.