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.
