Solving Mixed Integer Problems
This tutorial:
-
Shows how to solve mixed integer problems by enumeration and parallel optimization.
This method is supported by the fcmaes meta optimizer multiretry.
Some fcmaes optimizers also support discrete arguments directly, for example DE and MODE.
Others work well when continuous decision variables are mapped to integers using numpy’s
astype(int) or argsort(), for example BiteOpt.
Introduction
What if some variables in an optimization problem must be integers (see gear_train.py)? Or what if they must be integer multiples of some real value (see vessel.py)?
One option is to define all variables as continuous and map them to the nearest feasible value. That is what both examples above do. The vessel example also shows why derivative-based methods need care here. This mapping destroys smoothness in the objective function. All fcmaes algorithms are derivative-free, so this is usually not a problem.
Another option is to treat the integer variables as an integer sequence and enumerate all sequences. We can then solve the continuous optimization problem associated with each sequence.
Example: TandEM is formulated as a set of separate problems. Each one is associated with an integer sequence, namely the sequence of planets visited on the way to Titan and Enceladus. We could combine these problems into a single mixed integer formulation. Or we could optimize all sequences simultaneously. Let us compare the pros and cons of the three approaches for finding the best solution over all integer instances.
1) Solve all problems sequentially:
-
No additional variables means no extra complexity for each problem.
-
Sequential execution keeps memory requirements low.
-
Existing methods work out of the box.
But:
-
The number of sequences can be huge, which means high execution time.
-
We spend the same effort on bad sequences as on promising ones.
2) Mixed integer approach - add integer variables:
-
Memory requirements stay moderate, because we store data for only one optimization run.
-
It works well with derivative-free optimizers.
But:
-
The problem becomes much more complex because of the extra variables.
-
Integer sequences can only be restricted by boxed bounds. Sometimes domain knowledge lets us exclude additional sequences. For TandEM, boxed bounds give 24 possible sequences (
2*3*4), but ESA experts identified only 18 as feasible.
3) Simultaneous optimization of all sequences:
-
It can save execution time because bad sequences can be filtered out early.
-
Each individual problem keeps the same complexity as in the sequential approach.
But:
-
Data for all sequences must be stored in memory, which increases memory usage.
-
If there are many integer variables with wide bounds, the number of sequences can explode.
We conclude:
-
The sequential approach is best when all problem variants produce similar solutions.
-
The mixed integer approach is preferable when the number of sequences is very high and/or the underlying problem is not too complex.
-
The simultaneous optimization approach is preferable when domain knowledge can exclude sequences, the remaining number fits into memory, and/or the problem itself is very complex.
-
If the problem is very complex and the number of sequences is huge, we can start with the mixed integer approach to identify promising sequences quickly and then apply the simultaneous optimization approach to those sequences.
For TandEM we compare the mixed integer approach with the simultaneous optimization approach.
TandEm-MINLP
TandEM is an interplanetary mission problem that targets Titan and Enceladus, two moons of Saturn. The goal is to maximize the mass delivered to a high-eccentricity final orbit. We use the time-limited variant here because it is closer to a real mission scenario. With 18 parameters, TandEM is already hard to solve even when the planet sequence is fixed.
The departure from Earth and the arrival at Saturn are fixed, but there are three gravity-assist flybys in between. The original problem lists 18 suggested flyby sequences. We can reformulate it as TandEm-MINLP by adding three integer variables that describe the sequence. We restrict the three flybys to "Venus + Earth", "Venus + Earth + Mars" and "Venus + Earth + Mars + Jupiter" respectively. This yields 24 possible sequences and includes the 18 suggested ones.
TandEm is included in gtopx.cpp, and
astro.py provides its Tandem_minlp extension.
We apply coordinated retry in the usual way:
def test_optimizer(opt, problem, num_retries = 120000, num = 100,
value_limit = -10.0, log = logger()):
log.info(problem.name + ' ' + opt.name)
for _ in range(num):
ret = advretry.minimize(problem.fun, problem.bounds, value_limit,
num_retries, log, optimizer=opt)
test_optimizer(de_cma(1500), Tandem_minlp())and get the result described below.
|
Note
|
The original tandem-minlp.png result figure is not present in this repository.
|
TandEm-MINLP is tractable. We find the optimum at around 1500 km/s. Still, the complexity of this problem is comparable to Messenger-Full.
Eliminating sequences
Turning TandEm into a mixed integer problem is not the only option. On a 16-core CPU, this formulation already pushed coordinated retry close to its limits. Instead, we can optimize all sequence variants simultaneously. In each iteration, the worst 30% are removed. This shifts more computation toward the promising planet sequences.
def test_multiretry(num_retries = 512,
keep = 0.7, optimizer = de_cma(1500), logger = logger(), repeat = 50):
seqs = Tandem(0).seqs
n = len(seqs)
problems = [Tandem(i) for i in range(n)]
ids = [str(seqs[i]) for i in range(n)]
for _ in range(100):
# sort the problem variations
problem_stats = multiretry.minimize(problems, ids, num_retries, keep, optimizer, logger)
ps = problem_stats[0] # focus on the best one
for _ in range(repeat):
logger.info("problem " + ps.prob.name + ' ' + str(ps.id))
ps.retry(optimizer)
test_multiretry(repeat = 50)After 1326 sec we get the following sorted list of planet sequences:
[3, 2, 3, 3, 6] -1500.4667571840446 [3, 2, 2, 3, 6] -1328.751985466045 [3, 2, 2, 2, 6] -989.1379743688183 [3, 3, 4, 3, 6] -851.4661792233727 [3, 3, 2, 3, 6] -804.4525753263358 [3, 2, 3, 5, 6] -740.5926908771486 [3, 2, 4, 3, 6] -694.5843374074246 [3, 2, 2, 5, 6] -651.4881383543833 [3, 2, 3, 4, 6] -593.2141401623932 [3, 2, 4, 2, 6] -498.99456675910943 [3, 3, 4, 2, 6] -437.7147547067979 [3, 3, 4, 4, 6] -434.9099365281996 [3, 3, 3, 4, 6] -375.8823684957645 [3, 2, 2, 4, 6] -308.65460793463075 [3, 3, 2, 5, 6] -278.4845928699299 [3, 3, 4, 5, 6] -277.55722734268187 [3, 2, 4, 5, 6] -269.8725146307617 [3, 2, 3, 2, 6] -163.74018176257104 [3, 2, 4, 4, 6] -106.37445593961957 [3, 3, 2, 4, 6] -99.15011759793457 [3, 3, 2, 2, 6] -96.8058744309894 [3, 3, 3, 2, 6] -85.14366770071246 [3, 3, 3, 3, 6] -80.96971440986152 [3, 3, 3, 5, 6] -49.047704882833855
For the best sequence we already obtained an optimal solution. Otherwise we could continue by evaluating the top-ranked sequences further.
Cassini1-MINLP
At Midaco/GTOPX we find the example Cassini1-MINLP. Despite the name, this is not a linear programming problem. The objective function is nonlinear. MILP (Mixed Integer Linear Programming) problems have linear objective functions. They are easier to solve and require different algorithms for best efficiency.
It is recommended to read CEC2019 first, because it contains a detailed description of the problem.
We already mentioned the Cassini problem in Pykep gym results. Cassini1-MINLP extends that benchmark. The original Cassini Mission to Saturn used four gravity assists at Venus, Venus, Earth and Jupiter to save fuel. Cassini1-MINLP adds four variables that decide which planets are used for the four flybys. We show results for the problem from CEC2019 with the original integer variable bounds:
lb = [-1000.,30.,100.,30.,400.,1000., 1.0,1.0,1.0,1.0 ], ub = [0.,400.,470.,400.,2000.,6000., 9.0,9.0,9.0,9.0 ]
We also use limited bounds based on domain knowledge. For the first three flybys we allow only Mercury, Venus and Earth. For the last flyby we allow only Jupiter.
lb = [-1000.,30.,100.,30.,400.,1000., 0.51,0.51,0.51,4.51 ], ub = [0.,400.,470.,400.,2000.,6000., 3.49,3.49,3.49,5.49 ]
As we can see, both problems are tractable.
We also see that the best solution at around dv = 3.50 is protected by a dominant local optimum at dv = 3.63.
So another practical option is to try all sequences directly.
You cannot reproduce this result with fcmaes because there is a bug in the original result, and we reintroduced that bug here to make the comparison possible.
The issue is that Cassini1 has inequality constraints enforced by a penalty inside the GTOP code.
These constraints depend on the planet sequence, so the penalties also depend on the sequence.
CEC2019 unfortunately replaced them with fixed constraints that only work for the original planet sequence.
The resulting problem no longer matches the real mission.
Even so, it is still interesting because it creates a kind of worst-case scenario for the MINLP approach.
Cassini1-MINLP with correct constraints
Now let us apply fcmaes coordinated retry using the default optimizer de_cma
(see cassini_minlp.py).
We allow all nine planets in the sequence, which gives 9^4 = 6561 possible combinations:
from fcmaes.astro import Cassini1minlp
from fcmaes.optimizer import logger
from fcmaes.advretry import minimize
def test_optimizer(problem, num_retries = 100000, num = 100, value_limit = 100.0, log = logger()):
log.info(problem.name + ' ' + opt.name)
for i in range(num):
ret = minimize(problem.fun, problem.bounds, value_limit, num_retries, log)
test_optimizer(Cassini1minlp())We set num_retries = 100000 because we expect the problem to be hard.
The code was executed on an AMD 3950x 16-core processor at 4.0GHZ.
The results are surprising:
We observe strong resistance at deltaV = 2.896. After about 1000 sec, many runs improve on this value significantly. The best solution at deltaV = 2.17576 appears only once in 20 runs.
Solution 1: x = [-192.2248623732259, 47.13641665602527, 167.36749494279744, 315.1019785468309, 1095.7262611785677, 5695.214069302039, 3, 2, 3, 3] f(x) = 2.1757611862263295 Solution 2: x = [-738.4153386564027, 60.81046673324806, 162.5036882873382, 382.4428551935498, 1041.1379919643405, 1838.8750084727974, 3, 2, 2, 3] f(x) = 2.2509707061664836 Solution 3: x = [-749.7212652657448, 30.09190796959058, 175.81322285889598, 120.75196759589647, 1654.1828907793347, 1648.8003356144363, 3, 2, 1, 3] f(x) = 2.381303161294389
Fixing the objective function
As already mentioned in pykep gym results, there is a fundamental flaw in both the GTOP problems and their newer replacements pykep gym:
Restricting the coasting trajectory legs to single-revolution transfers has several nasty side effects, especially for the inner planets:
-
Global optima are much better shielded, because for long transfer times it is very unlikely to find a low-deltaV single-revolution transfer.
-
Many good solutions with multi-revolution coasting legs become invalid.
This flaw makes the GTOP problems artificially harder. It is not a real-world issue, because there is no reason to enforce the single-revolution restriction in practice.
Even if we do not know the best planet sequence in advance, we can still narrow the search space. For the first three flybys, going to any planet farther out than Earth makes little sense because it would slow the trajectory too much. By the same argument, the fourth encounter can be limited to Jupiter or a planet closer to the sun than Jupiter. This gives the following bounds, with equal continuous intervals assigned to each valid planet:
lb = [-1000.,30.,100.,30.,400.,1000., 0.51,0.51,0.51,2.51 ], ub = [0.,400.,470.,400.,2000.,6000., 3.49,3.49,3.49,5.49 ]
After implementing a fix for GTOP using the new PYKEP Lambert solver, we repeat the experiment with the restricted bounds that exclude the outer planets.
The curves are much smoother, which means there are many more good solutions.
The best one, with deltaV = 1.846 and planet sequence "Earth, Venus, Venus, Earth", is found in all ten runs after about 100 - 600 sec.
Why did they use Venus Venus Earth Jupiter for the real Cassini mission?
The choice of Venus, Venus, Earth, Jupiter for the real Cassini mission is easy to explain once we check the alternative planet sequences in the much more accurate model from pykep cassini2. In Pykep gym results we show a good solution for the original planet sequence with deltaV = 0.729 km/s. A quick 5 min check for the other sequences gives:
-
VVEJ: deltaV = 0.729 km/s
-
EVEJ: deltaV = 2.434 km/s
-
EEVE: deltaV = 2.606 km/s
-
EVVE: deltaV = 3.534 km/s
-
EMMJ: deltaV = 9.571 km/s
Using deep-space maneuvers and a more accurate model of the real mission shows that VVEJ needs only a fraction of the deltaV required by the other planet sequences.
GTOC1
The GTOC1 competition was already discussed in detail in the Tutorials.
Here we focus on finding a good planet sequence.
We use the "fixed multirevolution Lambert" version of the GTOP code, which is not yet included in fcmaes.
For GTOC1 this is essential, because it enables solutions that improve on the winning one from JPL.
We also need a long planet sequence, so we added a time of flight ⇐ 30 year constraint to the code.
Our boxed constraints for the timing variables are:
lb = [3000.,14.,14.,14.,14.,14.,14.,300.,300.,300.,300.],
ub = [10000.,2000.,2000.,2000.,2000.,2000.,2000.,3000.,1000.,4000.,1000.]
A typical planet sequence is [3, 2, 2, 3, 2, 2, 3, 3, 6, 5, 10], with ten trajectory legs starting at Earth (3) and ending at the asteroid that threatens our planet (10).
Using domain knowledge, we fix the last two flybys to Saturn and Jupiter and allow only Venus and Earth flybys in between.
That leaves 2^7 = 128 different planet sequences.
We do not expect large differences between sequences, so instead of using MINLP we choose the simultaneous optimization approach:
def sequences():
for p1 in range(2,4):
for p2 in range(2, 4):
for p3 in range(2, 4):
for p4 in range(2, 4):
for p5 in range(2, 4):
for p6 in range(2, 4):
for p7 in range(2, 4):
yield[p1,p2,p3,p4,p5,p6,p7]
def test_multiretry(num_retries = 1024),
keep = 0.7, optimizer = de_cma(1500), logger = logger(), repeat = 100):
problems = []
ids = []
for seq in sequences():
problems.append(Gtoc1(planets = seq))
ids.append(str(seq))
for _ in range(100):
problem_stats = multiretry.minimize(problems, ids, num_retries, keep, optimizer, logger)
ps = problem_stats[0]
for _ in range(repeat):
logger.info("problem " + ps.prob.name + ' ' + str(ps.id))
ps.retry(optimizer)Because of the enormous complexity of this task, we use a base retry count of 1024 for each iteration. In every iteration, 30% of the problem variants and sequences are filtered out. After about 25 minutes on a 3950X 16-core CPU, we get the first ranking:
[3, 3, 3, 3, 2, 3, 3] -1708989.6693652852 [3, 2, 3, 2, 3, 3, 3] -1696396.0489073088 [3, 3, 3, 2, 3, 3, 3] -1682620.699099308 [3, 3, 2, 3, 3, 3, 3] -1680375.8026656334 [2, 2, 3, 2, 2, 3, 3] -1663756.207142463 [3, 2, 3, 3, 2, 3, 3] -1659601.1495679119 [3, 2, 3, 2, 2, 3, 3] -1655046.4219106187 [3, 3, 3, 2, 2, 3, 3] -1639178.8846807072 [3, 3, 2, 2, 3, 2, 3] -1637261.6500321375 [3, 3, 2, 2, 3, 3, 3] -1636352.5885824538 ...
But the final result of the multiretry sorting process appears only after 8 more hours:
[2, 2, 3, 2, 3, 3, 3] -1847110.3065493396 ----- [3, 2, 3, 2, 3, 3, 3] -1847069.2810261745 ----- [2, 3, 2, 3, 2, 3, 3] -1843218.757455708 ----- [3, 2, 3, 3, 2, 3, 3] -1836566.1917796042 [2, 3, 2, 2, 3, 3, 3] -1827985.9268557034 [3, 2, 3, 2, 2, 3, 3] -1816419.265119756 [3, 3, 2, 2, 3, 3, 3] -1813416.482933015 [3, 3, 2, 3, 2, 3, 3] -1813359.2309791856 [3, 3, 3, 2, 3, 3, 3] -1782243.9263929802 [3, 2, 3, 3, 3, 3, 3] -1780021.7413547586 [3, 3, 2, 2, 2, 3, 3] -1779001.0378138914 [2, 3, 2, 2, 2, 3, 3] -1770922.5899577655 [3, 2, 2, 3, 2, 3, 3] -1765913.2291847111 [3, 2, 2, 2, 3, 3, 3] -1760166.55383832 [2, 2, 3, 2, 2, 3, 3] -1760078.2693866405 [3, 3, 3, 3, 2, 3, 3] -1759074.1787529092 [2, 3, 3, 3, 2, 3, 3] -1740435.1296423872 [3, 3, 2, 2, 3, 2, 3] -1733650.8112690973 [3, 3, 2, 3, 3, 3, 3] -1727629.2676451143 [3, 3, 3, 2, 2, 3, 3] -1712522.4017491455 [2, 3, 3, 2, 3, 3, 3] -1711634.8412391657 [3, 2, 3, 3, 2, 2, 3] -1710975.450628684 ...
Further evaluation of the first three sequences with regular coordinated retry gives:
[3,2,3,2,2,3,3,3,6,5,10] score = -1927767 [3,2,3,2,3,2,3,3,6,5,10] score = -1924321 [3,3,2,3,2,3,3,3,6,5,10] score = -1906362
The simultaneous retry predicted the order of the first three sequences correctly. Most likely, any of them would have been good enough to win the competition, since JPL scored only 1850000. All three have very low deltaV impulses, so converting them into a low-thrust trajectory is not difficult.