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. ]) 1 Answer
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.
3
