11: Optimalizační problémy
Tento notebook je výukovým materiálem v předmětu BI-JUL.21 vyučovaném v zimním semestru akademického roku 2021/2022 Tomášem Kalvodou. Tvorba těchto materiálů byla podpořena NVS FIT.
Hlavní stránkou předmětu, kde jsou i další notebooky a zajímavé informace, je jeho Course Pages stránka.
1. Úvod
Nejprve je nutné učinit několik úvodních vysvětlujících poznámek. Co přesně znamená slovní spojení "optimalizačním problémem"?
1.1 Matematická formulace optimalizačního problému
Pod "optimalizačním problémem" máme na mysli následující matematickou úlohu:
Najděte extrém (minimum či maximum) funkce , kde množina je popsána pomocí funkcí , , a , , následovně
Mluvíme o "rovnostních" a "nerovnostních" podmínkách, v některých problémech mohou být podmínky jen jednoho druhu nebo dokonce nemusíme mít podmínky žádné.
Množina může být otevřená podmnožina , nebo třeba podmnožina , či dokonce podmnožina .
Množina je množina tzv. "přípustných" řešení. Součástí úlohy je i nalezení nějakého jejího prvku (je neprázdná?). Přípustné řešení splňující (případně ) se nazývá "optimální řešení" a "optimální hodnota" daného problému (optimální hodnota je právě jedna - existuje-li; optimálních řešení může být více).
Nejčastější kategorizace těchto problémů je podle charakteru používaných funkcí a charakteru proměnných, používají se tyto zkratky:
- MI: Mixed Integer, některé proměnné mohou být omezeny na celočíselné hodnoty (pouze binární hodnoty 0 a 1 se většionu počítají do této kategorie),
- LP: Linear programming, objektivní funkce a omezení jsou lineární (affinní) funkce,
- QP: Quadratic programming, objektivní funkce kvadratická, ale omezení stále lineární (affinní) funkce,
- CP: Convex programming, objektivní funkce a nerovnostní omezení konvexní, rovnostní omezení lineární (affinní) funkce.
1.2 Historická perspektiva
Optimalizační úlohy se nejprve objevily ve fyzice. Příroda má totiž tendenci minimalizovat různé veličiny (energii, čas, vzdálenost,...). Z computer science perspektivy je první zajímavý SIMPLEX algoritmus (Dantzig, 40. léta dvacátého století) řešící optimalizační úlohu s lineárními objektivními funkcemi a omezeními. SIMPLEX současně dává algoritmické řešení soustavy lineárních nerovnic. Optimalizační úlohy během dvacátého století vzkvétaly a nacházely další uplatnění například ve strojovém učení a plánování (v tomto smyslu se zde používá slovíčko programming). Často už ovšem nejde pouze o lineární problémy!
1.3 Julia nástroje
V Julia je asi nejpoužívanějším nástrojem na řešení optimalizačních úloh balíček JuMP.jl
.
Jde o zkratku Julia for Mathematical Programming.
Tento balíček poskytuje výrazivo umožňující optimalizační problém snadno formulovat.
K samotnému řešení problému pak používá některou z existujících knihoven nebo Julia balíčků.
Těch přirozeně existuje celá řada a každá umí řešit různé typy problémů, JuMP se je snaží všechny zastřešit a odstínit uživatele od práce s API samotných knihoven/balíčků.
Některé z těchto knihoven jsou open source, některé jsou komerční.
Kompletní přehled podporovaných "řešičů" naleznete v dokumentaci zde.
K zprovoznění tohoto notebooku budeme potřebovat balíček JuMP.jl
samotný, nainstalujte ho standardně pomocí ] add JuMP
a poté nainstalujte i některý z "řešičů" (viz níže).
using JuMP, PyPlot
JuMP.jl
2. Základní formulace problému ("modelu") v Formulace problému (také bývá zvykem hovořit o "modelu") v JuMP velmi věrně kopíruje jeho matematickou strukturu. Musíme vytvořit proměnné, omezení, objektivní funkci a případně nastavit další parametry. JuMP k tomu velmi šikovně využívá Julia makra. Poté, co takovýto model vytvoříme, může ho JuMP samostatně poslat zvolenému řešiči. O tuto práci se už jako uživatelé starat nemusíme a JuMP ji udělá za nás (není to triviální, každý řešič může mít jiné API, JuMP se snaží i naše funkce optimalizovat nebo i dopočítává derivace, pokud je řešič vyžaduje).
Ukažme si tento proces na jednoduchém a snadno představitelném příkladě. Zájemcům doporučuji projít si velmi dobře sepsanou dokumentaci JuMP.
Minimalizujme funkci na obdélníku . Tuto situaci si můžeme velmi snadno vizualizovat (vzhledem k malému počtu proměnných).
xs = LinRange(-1, 1, 50)
ys = LinRange(-1, 1, 50)
fs = [ x^2 - x*y - y^2 for y in ys, x in xs ]
plt.xlabel("x")
plt.ylabel("y")
plot_surface(xs, ys, fs);
Ačkoliv tento graf vypadá atraktivně, není moc přehledný. Daleko lepší je contour plot nebo heat map.
plt.xlabel("x")
plt.ylabel("y")
contourf(xs, ys, fs, 20)
plt.colorbar()
PyObject <matplotlib.colorbar.Colorbar object at 0x7f09cc092e80>
Z grafů se zdá, že minimum bude nabyto někde na hranici obdélníku. Pojďme to prověřit.
Nejprve musíme vytvořit Model
, při tom můžeme specifikovat jaký řešič se použije.
Nejprve ale musíme importovat balíček s řešičem, zde Ipopt, (] add Ipopt
).
using Ipopt
Poté můžeme vytvořit náš model a předat mu řešič.
model = Model(Ipopt.Optimizer);
Náš problém má dvě proměnné a .
Ukážeme si dva způsoby jak definovat omezení.
K definici proměnných se používá makro @variable
, kterému můžeme přidat i jednoduchá omezení na rozsah:
@variable(model, -1 <= x <= 1)
@variable(model, y);
Případná další, potenciální výrazně komplikovanější, omezení můžeme definovat pomocí makra @constraint
.
Omezme pomocí něho tedy i druhou proměnnou (proměnná by v těchto zápisech měla být na levé straně nerovnosti).
@constraint(model, y >= -1)
@constraint(model, y <= 1);
V proměnné x
je uložen JuMP objekt typu VariableRef
, který nám umožňuje vytvářet různé symbolické výrazy a pomocí něhož poté budeme získávat řešení.
typeof(x)
VariableRef
2 * x
A konečně k definici objektivní funkce slouží makro @objective
:
@objective(model, Min, x^2 - x*y - y^2);
Makro @objective
si poradí s lineární nebo kvadratickou objektivní funkcí. Pokud bychom chtěli pracovat s "výrazněji" nelineárními podmínkami nebo objektivními funkcemi, musíme explicitně použít makra @NLconstraint
a @NLobjective
.
Než se pustíme do řešení, můžeme si nechat model zobrazit (toto se hodí jen pro menší modely):
print(model)
latex_formulation(model)
Případně:
model
A JuMP Model Minimization problem with: Variables: 2 Objective function type: QuadExpr `AffExpr`-in-`MathOptInterface.GreaterThan{Float64}`: 1 constraint `AffExpr`-in-`MathOptInterface.LessThan{Float64}`: 1 constraint `VariableRef`-in-`MathOptInterface.GreaterThan{Float64}`: 1 constraint `VariableRef`-in-`MathOptInterface.LessThan{Float64}`: 1 constraint Model mode: AUTOMATIC CachingOptimizer state: EMPTY_OPTIMIZER Solver name: Ipopt Names registered in the model: x, y
Řešení spustíme metodou optimize!
:
optimize!(model)
This is Ipopt version 3.13.4, running with linear solver mumps. NOTE: Other linear solvers might be more efficient (see Ipopt documentation). Number of nonzeros in equality constraint Jacobian...: 0 Number of nonzeros in inequality constraint Jacobian.: 2 Number of nonzeros in Lagrangian Hessian.............: 3 Total number of variables............................: 2 variables with only lower bounds: 0 variables with lower and upper bounds: 1 variables with only upper bounds: 0 Total number of equality constraints.................: 0 Total number of inequality constraints...............: 2 inequality constraints with only lower bounds: 1 inequality constraints with lower and upper bounds: 0 inequality constraints with only upper bounds: 1 iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls 0 0.0000000e+00 0.00e+00 0.00e+00 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0 1 -2.3261961e-34 0.00e+00 2.00e-07 -1.7 1.42e-17 0.0 1.00e+00 1.00e+00 0 2 -4.7796348e-34 0.00e+00 1.50e-09 -3.8 5.79e-18 0.4 1.00e+00 1.00e+00 0 3 -3.6032081e-32 0.00e+00 1.84e-11 -5.7 1.51e-16 -0.1 1.00e+00 1.00e+00 0 4 -8.0681656e-32 0.00e+00 2.53e-14 -8.6 8.49e-17 0.4 1.00e+00 1.00e+00T 0 Number of Iterations....: 4 (scaled) (unscaled) Objective...............: -8.0681656006255618e-32 -8.0681656006255618e-32 Dual infeasibility......: 2.5260350807739885e-14 2.5260350807739885e-14 Constraint violation....: 0.0000000000000000e+00 0.0000000000000000e+00 Complementarity.........: 2.5059035596802921e-09 2.5059035596802921e-09 Overall NLP error.......: 2.5059035596802921e-09 2.5059035596802921e-09 Number of objective function evaluations = 5 Number of objective gradient evaluations = 5 Number of equality constraint evaluations = 0 Number of inequality constraint evaluations = 5 Number of equality constraint Jacobian evaluations = 0 Number of inequality constraint Jacobian evaluations = 1 Number of Lagrangian Hessian evaluations = 1 Total CPU secs in IPOPT (w/o function evaluations) = 0.003 Total CPU secs in NLP function evaluations = 0.000 EXIT: Optimal Solution Found.
Tento log výpočtu nepochází z JuMP, ale přímo z Ipopt, nebo i jiných řešičů. Potlačit ho můžeme zavolání set_silent(model)
.
Pojďme prozkoumat výsledek:
solution_summary(model)
* Solver : Ipopt * Status Termination status : LOCALLY_SOLVED Primal status : FEASIBLE_POINT Dual status : FEASIBLE_POINT Message from the solver: "Solve_Succeeded" * Candidate solution Objective value : -8.068165600625562e-32 * Work counters Solve time (sec) : 0.00427
K optimálnímu řešení se dostaneme pomocí metody value
.
value(x)
9.167787467101626e-17
value(y)
2.5613406954838343e-16
value.([x, y])
2-element Vector{Float64}: 9.167787467101626e-17 2.5613406954838343e-16
A k optimální hodnotě pomocí objective_value
:
objective_value(model)
-8.068165600625562e-32
Vidíme, že se Ipopt řešič nedokázal dostat ze stacionárního bodu , kde pravděpodobně začal. Tento bod dokonce není ani bodem lokálního extrému! Musíme do řešiče malinko kopnout (alternativně by mohlo pomoci změnit optimalizační metodu / zkusit jiný řešič).
set_start_value(x, 0.23)
Optimalizujme znovu!
optimize!(model)
This is Ipopt version 3.13.4, running with linear solver mumps. NOTE: Other linear solvers might be more efficient (see Ipopt documentation). Number of nonzeros in equality constraint Jacobian...: 0 Number of nonzeros in inequality constraint Jacobian.: 2 Number of nonzeros in Lagrangian Hessian.............: 3 Total number of variables............................: 2 variables with only lower bounds: 0 variables with lower and upper bounds: 1 variables with only upper bounds: 0 Total number of equality constraints.................: 0 Total number of inequality constraints...............: 2 inequality constraints with only lower bounds: 1 inequality constraints with lower and upper bounds: 0 inequality constraints with only upper bounds: 1 iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls 0 5.2900000e-02 0.00e+00 4.60e-01 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0 1 1.0828786e-02 0.00e+00 9.04e-02 -1.0 9.04e-02 0.0 1.00e+00 1.00e+00f 1 2 -3.9823732e-03 0.00e+00 1.12e-01 -2.5 4.21e-02 0.4 1.00e+00 1.00e+00f 1 3 -6.2968265e-01 0.00e+00 5.66e-01 -2.5 6.36e-01 -0.1 1.00e+00 1.00e+00f 1 4 -1.2314300e+00 0.00e+00 7.84e-01 -2.5 3.52e-01 0.4 1.00e+00 8.11e-01f 1 5 -1.2289751e+00 0.00e+00 1.54e+00 -2.5 1.05e-01 - 1.00e+00 3.12e-02f 6 6 -1.2472887e+00 0.00e+00 2.83e-08 -2.5 1.23e-01 - 1.00e+00 1.00e+00f 1 7 -1.2498526e+00 0.00e+00 1.50e-09 -3.8 2.23e-03 - 1.00e+00 1.00e+00f 1 8 -1.2499982e+00 0.00e+00 1.84e-11 -5.7 1.41e-04 - 1.00e+00 1.00e+00f 1 9 -1.2500000e+00 0.00e+00 2.53e-14 -8.6 1.64e-06 - 1.00e+00 1.00e+00f 1 Number of Iterations....: 9 (scaled) (unscaled) Objective...............: -1.2500000224957013e+00 -1.2500000224957013e+00 Dual infeasibility......: 2.5313084961453569e-14 2.5313084961453569e-14 Constraint violation....: 0.0000000000000000e+00 0.0000000000000000e+00 Complementarity.........: 2.5121111717748235e-09 2.5121111717748235e-09 Overall NLP error.......: 2.5121111717748235e-09 2.5121111717748235e-09 Number of objective function evaluations = 19 Number of objective gradient evaluations = 10 Number of equality constraint evaluations = 0 Number of inequality constraint evaluations = 19 Number of equality constraint Jacobian evaluations = 0 Number of inequality constraint Jacobian evaluations = 1 Number of Lagrangian Hessian evaluations = 1 Total CPU secs in IPOPT (w/o function evaluations) = 0.006 Total CPU secs in NLP function evaluations = 0.000 EXIT: Optimal Solution Found.
Optimální řešení:
value.([x, y])
2-element Vector{Float64}: 0.500000002821662 1.0000000089982806
Kdybychom se vychýlili na druhou stranu, tak dokonvergujeme do jiného lokálního minima:
set_start_value(x, -0.1)
optimize!(model)
This is Ipopt version 3.13.4, running with linear solver mumps. NOTE: Other linear solvers might be more efficient (see Ipopt documentation). Number of nonzeros in equality constraint Jacobian...: 0 Number of nonzeros in inequality constraint Jacobian.: 2 Number of nonzeros in Lagrangian Hessian.............: 3 Total number of variables............................: 2 variables with only lower bounds: 0 variables with lower and upper bounds: 1 variables with only upper bounds: 0 Total number of equality constraints.................: 0 Total number of inequality constraints...............: 2 inequality constraints with only lower bounds: 1 inequality constraints with lower and upper bounds: 0 inequality constraints with only upper bounds: 1 iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls 0 1.0000000e-02 0.00e+00 2.00e-01 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0 1 2.0086300e-03 0.00e+00 3.99e-02 -1.0 3.99e-02 0.0 1.00e+00 1.00e+00f 1 2 -7.3425290e-04 0.00e+00 4.79e-02 -2.5 1.80e-02 0.4 1.00e+00 1.00e+00f 1 3 -1.1931993e-01 0.00e+00 2.46e-01 -2.5 2.77e-01 -0.1 1.00e+00 1.00e+00f 1 4 -2.6651816e-01 0.00e+00 3.66e-01 -2.5 1.54e-01 0.4 1.00e+00 1.00e+00f 1 5 -1.2275584e+00 0.00e+00 7.84e-01 -2.5 4.86e+01 -0.1 3.01e-02 1.10e-02f 1 6 -1.2298807e+00 0.00e+00 1.84e+00 -2.5 1.66e-01 - 1.00e+00 1.20e-02f 2 7 -1.2472209e+00 0.00e+00 2.83e-08 -2.5 1.33e-01 - 1.00e+00 1.00e+00f 1 8 -1.2498528e+00 0.00e+00 1.50e-09 -3.8 2.31e-03 - 1.00e+00 1.00e+00f 1 9 -1.2499982e+00 0.00e+00 1.84e-11 -5.7 1.42e-04 - 1.00e+00 1.00e+00f 1 iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls 10 -1.2500000e+00 0.00e+00 2.51e-14 -8.6 1.64e-06 - 1.00e+00 1.00e+00f 1 Number of Iterations....: 10 (scaled) (unscaled) Objective...............: -1.2500000224957013e+00 -1.2500000224957013e+00 Dual infeasibility......: 2.5059035640133008e-14 2.5059035640133008e-14 Constraint violation....: 0.0000000000000000e+00 0.0000000000000000e+00 Complementarity.........: 2.5121141971489301e-09 2.5121141971489301e-09 Overall NLP error.......: 2.5121141971489301e-09 2.5121141971489301e-09 Number of objective function evaluations = 13 Number of objective gradient evaluations = 11 Number of equality constraint evaluations = 0 Number of inequality constraint evaluations = 13 Number of equality constraint Jacobian evaluations = 0 Number of inequality constraint Jacobian evaluations = 1 Number of Lagrangian Hessian evaluations = 1 Total CPU secs in IPOPT (w/o function evaluations) = 0.007 Total CPU secs in NLP function evaluations = 0.000 EXIT: Optimal Solution Found.
Optimální řešení:
value.([x,y])
2-element Vector{Float64}: -0.5000000028216589 -1.0000000089982806
objective_value(model)
Najít globální minima a maxima je těžká úloha, pokud o daných funkcích nic více nevíme (typicky potřebujeme alespoň konvexitu). Řešiče v jednom běhu hledají lokální extrémy.
Pro demonstraci ukažme použití i komerčního nástroje Gurobi (pro akademické účely volně k použití; je potřeba nainstalovat balíček ] add Gurobi
a získat licenci, viz oficiální web):
using Gurobi
model = Model(Gurobi.Optimizer)
set_optimizer_attribute(model, "NonConvex", 2)
@variable(model, -1 <= x <= 1)
@variable(model, -1 <= y <= 1)
@objective(model, Min, x^2 - x*y - y^2);
optimize!(model)
Set parameter Username Academic license - for non-commercial use only - expires 2022-11-07 Set parameter NonConvex to value 2 Set parameter NonConvex to value 2 Gurobi Optimizer version 9.5.0 build v9.5.0rc5 (linux64) Thread count: 4 physical cores, 8 logical processors, using up to 8 threads Optimize a model with 0 rows, 2 columns and 0 nonzeros Model fingerprint: 0x29307fce Model has 3 quadratic objective terms Coefficient statistics: Matrix range [0e+00, 0e+00] Objective range [0e+00, 0e+00] QObjective range [2e+00, 2e+00] Bounds range [1e+00, 1e+00] RHS range [0e+00, 0e+00] Continuous model is non-convex -- solving as a MIP Found heuristic solution: objective -1.0000000 Presolve time: 0.00s Presolved: 0 rows, 2 columns, 0 nonzeros Presolved model has 3 quadratic objective terms Variable types: 1 continuous, 1 integer (1 binary) Root relaxation: objective -1.250000e+00, 2 iterations, 0.00 seconds (0.00 work units) Nodes | Current Node | Objective Bounds | Work Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time * 0 0 0 -1.2500000 -1.25000 0.00% - 0s Explored 1 nodes (2 simplex iterations) in 0.00 seconds (0.00 work units) Thread count was 8 (of 8 available processors) Solution count 1: -1.25 No other solutions better than -1.25 Optimal solution found (tolerance 1.00e-04) Best objective -1.250000000000e+00, best bound -1.250000000000e+00, gap 0.0000% User-callback calls 91, time in user-callback 0.00 sec
Opět získáváme jedno z možných optimálních řešení:
value.([x,y])
2-element Vector{Float64}: -0.5 -1.0
objective_value(model)
-1.25
3. Cvičení / Úlohy
Po tomto rychlém úvodu do JuMP si ukážeme několik zajímavých a zejména různorodých optimalizačních úloh.
3.1 "Tradiční" zemědělská úloha
Takřka v každé učebnici lineárního programování narazíte na úlohu následujícího typu:
Mějme farmáře, který obdělává pole o rozloze . Pole plánuje osít třemi plodinami (1., 2. a 3., nebudeme se pouštět do konkrétního pojmenování).
- Cena osetí jednotky plochy danou plodinou je popořadě , a .
- Péče o jednotku plochy oseté danou plodinou ho bude stát , a .
- Výnos z jednotky plochy oseté danou plodinou bude , a .
Jak má náš zemědělec využít své pole, aby maximalizoval zisk? Ani v tomto kurzu si nedovolíme tento typ úlohy nezmínit!
Rozeberme situaci podrobně. Neznámými jsou zřejmě plochy oseté jednotlivými plodinami. Označme si je , a . Jistě musí platit následující čtyři nerovnosti:
Snažíme se maximalizovat objektivní funkci (profit; k zjednodušení zápisu jsme přešli k vektorové notaci)
Skutečně jde o lineární program: omezující funkce i objektivní funkce jsou lineární funkce (stačí affinní -- tj. lineární plus konstanta).
Poznámka: Úlohu by šlo samozřejmě dále komplikovat, například přidat horní limit rozpočtu atp. Analogické úlohy lze formulovat také pro dietologické problémy (nastavení jídelníčku s požadovými nutričními požadavky). Opět vidíme, v že název lineární programování by bylo bývalo možná lepší překládat jako lineární plánování.
Zafixujme nějaké konkrétní hodnoty a pokusme se problém vyřešit pomocí Ipopt a JuMP.
using JuMP, Ipopt
S = 10
v1, v2, v3 = 5.5, 6, 6.5
p1, p2, p3 = 1, 1.5, 2
o1, o2, o3 = 1, 1, 1.5;
Sestavení modelu:
model = Model(Ipopt.Optimizer)
@variable(model, x[1:3] >= 0)
@constraint(model, x[1] + x[2] + x[3] <= 10)
@objective(model, Max, (v1 - p1 - o1) * x[1] + (v2 - p2 - o2) * x[2] + (v3 - p3 - o3) * x[3])
print(model)
Spusťme řešení!
optimize!(model)
This is Ipopt version 3.13.4, running with linear solver mumps. NOTE: Other linear solvers might be more efficient (see Ipopt documentation). Number of nonzeros in equality constraint Jacobian...: 0 Number of nonzeros in inequality constraint Jacobian.: 3 Number of nonzeros in Lagrangian Hessian.............: 0 Total number of variables............................: 3 variables with only lower bounds: 3 variables with lower and upper bounds: 0 variables with only upper bounds: 0 Total number of equality constraints.................: 0 Total number of inequality constraints...............: 1 inequality constraints with only lower bounds: 0 inequality constraints with lower and upper bounds: 0 inequality constraints with only upper bounds: 1 iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls 0 9.9999900e-02 0.00e+00 2.50e+00 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0 1 1.4299999e+00 0.00e+00 3.46e+00 -1.0 3.99e-01 - 2.22e-01 1.00e+00f 1 2 3.4483052e+01 0.00e+00 2.68e+00 -1.0 3.60e+01 - 1.85e-02 2.63e-01f 1 3 3.4711502e+01 0.00e+00 2.25e+00 -1.0 7.56e-01 - 1.00e+00 1.66e-01f 2 4 3.4804478e+01 0.00e+00 1.00e-06 -1.0 4.65e-02 - 1.00e+00 1.00e+00f 1 5 3.4990149e+01 0.00e+00 2.83e-08 -2.5 1.79e-01 - 1.00e+00 1.00e+00f 1 6 3.4999698e+01 0.00e+00 1.50e-09 -3.8 1.26e-02 - 1.00e+00 1.00e+00f 1 7 3.4999997e+01 0.00e+00 1.84e-11 -5.7 3.00e-04 - 1.00e+00 1.00e+00f 1 8 3.5000000e+01 0.00e+00 2.53e-14 -8.6 3.70e-06 - 1.00e+00 1.00e+00f 1 Number of Iterations....: 8 (scaled) (unscaled) Objective...............: -3.5000000349986635e+01 3.5000000349986635e+01 Dual infeasibility......: 2.5313084961453569e-14 2.5313084961453569e-14 Constraint violation....: 0.0000000000000000e+00 0.0000000000000000e+00 Complementarity.........: 2.5072641342280509e-09 -2.5072641342280509e-09 Overall NLP error.......: 2.5072641342280509e-09 2.5313084961453569e-14 Number of objective function evaluations = 11 Number of objective gradient evaluations = 9 Number of equality constraint evaluations = 0 Number of inequality constraint evaluations = 11 Number of equality constraint Jacobian evaluations = 0 Number of inequality constraint Jacobian evaluations = 1 Number of Lagrangian Hessian evaluations = 1 Total CPU secs in IPOPT (w/o function evaluations) = 0.005 Total CPU secs in NLP function evaluations = 0.000 EXIT: Optimal Solution Found.
Optimální řešení a hodnota:
value.(x)
3-element Vector{Float64}: 5.000000052134667 5.0000000521347765 0.0
objective_value(model)
35.000000349986635
Naše doporučení zemědělci by bylo: osej pomocí první a druhé plodiny půl na půl. Očekávaný zisk pak budeš mít "".
Když se podíváte na vstupní data, tak toto zcela jistě není na první pohled patrné!
3.2 Řetězovka
Představte si řetěz délky , který je zavěšen v bodech a , , v homogenním gravitačním poli (v záporném směru osy ). Jak se tento řetěz prověsí? Představte si třeba různé řetězové ploty nebo závěsy mostů. Fyzika (klasická mechanika, detaily vynecháváme) bude tlačit řetěz do stavu s minimální potenciální energií. Pojďme tuto minimalizaci popsat pomocí JuMP.
using JuMP, Ipopt, Gurobi
Řetěz budeme idealizovat jako hmotných bodů spojených nehmotným segmentem. Zajímá nás samozřejmě limita , resp. její aproximace. Zvolíme následující parametry:
n = 10
D = 5
a = 2
L = 8
Δ = L / n
if D > L
@warn "Tohle neprověsíš!"
end
Poté definujeme model. Všimněte si, že objektivní funkce je lineární a netriviální podmínky jsou kvadratické, jde ale o nelineární problém.
model = Model(Ipopt.Optimizer)
# set_optimizer_attribute(model, "NonConvex", 2) # Gurobi
# neznámé jsou souřadnice bodů (x[i], y[i]), i = 1,...,n
@variable(model, x[1:n])
@variable(model, y[1:n])
# okrajové podmínky
@constraints(model, begin
x[1] == 0
x[n] == D
y[1] == a
y[n] == a
end)
# spojky řetězu, resp. vzdálenosti sousedních bodů
for i=1:(n-1)
@constraint(model, (x[i] - x[i+1])^2 + (y[i] - y[i+1])^2 == Δ^2)
end
# objektivní funkce: vynecháváme fyzikální konstanty, které by jen škálovaly celkovou
# funkční hodnotu (krát hmotnost hmotného bodu krát gravitační zrychlení).
# Efekt těžšího dílu řetězu by se projevil násobením příslušeného členu větší konstantou než 1.
@objective(model, Min, sum(y[i] for i=1:n));
print(model)
Pokusme se nalézt řešení.
optimize!(model)
****************************************************************************** This program contains Ipopt, a library for large-scale nonlinear optimization. Ipopt is released as open source code under the Eclipse Public License (EPL). For more information visit https://github.com/coin-or/Ipopt ****************************************************************************** This is Ipopt version 3.13.4, running with linear solver mumps. NOTE: Other linear solvers might be more efficient (see Ipopt documentation). Number of nonzeros in equality constraint Jacobian...: 76 Number of nonzeros in inequality constraint Jacobian.: 0 Number of nonzeros in Lagrangian Hessian.............: 54 Total number of variables............................: 20 variables with only lower bounds: 0 variables with lower and upper bounds: 0 variables with only upper bounds: 0 Total number of equality constraints.................: 13 Total number of inequality constraints...............: 0 inequality constraints with only lower bounds: 0 inequality constraints with lower and upper bounds: 0 inequality constraints with only upper bounds: 0 iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls 0 0.0000000e+00 5.00e+00 1.00e+00 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0 1 -7.9600000e+02 1.04e+04 2.32e+10 -1.0 1.00e+02 -2.0 1.00e+00 1.00e+00f 1 2 -7.9474186e+02 2.61e+03 9.52e+10 -1.0 2.62e+01 9.8 1.00e+00 1.00e+00h 1 3 -7.8410906e+02 6.53e+02 2.06e+11 -1.0 2.02e+01 9.4 1.00e+00 1.00e+00h 1 4 -7.8320856e+02 1.95e+02 3.09e+11 -1.0 1.13e+01 10.7 1.00e+00 1.00e+00h 1 5 -7.7907489e+02 6.67e+01 1.57e+11 -1.0 1.05e+01 10.2 1.00e+00 1.00e+00h 1 6 -7.6201336e+02 7.24e+01 5.34e+10 -1.0 8.76e+00 9.7 1.00e+00 1.00e+00h 1 7 -6.9753737e+02 6.88e+01 1.91e+10 -1.0 9.62e+00 9.3 1.00e+00 1.00e+00h 1 8 -4.7690385e+02 4.69e+01 1.38e+10 -1.0 2.25e+01 8.8 1.00e+00 1.00e+00h 1 9 1.9675139e+01 1.56e+00 3.30e+09 -1.0 4.32e+02 - 1.00e+00 1.00e+00H 1 iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls 10 1.9988291e+01 1.43e+00 5.89e+07 -1.0 1.36e-01 - 1.00e+00 1.00e+00h 1 11 1.9999923e+01 1.45e+00 3.56e+07 -1.0 6.21e-02 - 1.00e+00 1.00e+00h 1 12 2.0000026e+01 1.42e+00 6.04e+06 -1.0 2.80e-02 - 1.00e+00 1.00e+00h 1 13 2.0000026e+01 1.39e+00 3.94e+06 -1.0 3.05e-02 - 1.00e+00 1.00e+00h 1 14 2.0000025e+01 1.36e+00 2.51e+06 -1.0 3.48e-02 - 1.00e+00 1.00e+00h 1 15 2.0000024e+01 1.32e+00 2.48e+06 -1.0 3.68e-02 - 1.00e+00 1.00e+00h 1 16 2.0000024e+01 1.28e+00 2.44e+06 -1.0 3.80e-02 - 1.00e+00 1.00e+00h 1 17 2.0000023e+01 1.25e+00 2.41e+06 -1.0 3.85e-02 - 1.00e+00 1.00e+00h 1 18 2.0000022e+01 1.21e+00 2.37e+06 -1.0 3.87e-02 - 1.00e+00 1.00e+00h 1 19 2.0000021e+01 1.17e+00 2.33e+06 -1.0 3.87e-02 - 1.00e+00 1.00e+00h 1 iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls 20 2.0000021e+01 1.13e+00 2.29e+06 -1.0 3.85e-02 - 1.00e+00 1.00e+00h 1 21 2.0000020e+01 1.09e+00 2.25e+06 -1.0 3.83e-02 - 1.00e+00 1.00e+00h 1 22 2.0000019e+01 1.05e+00 2.20e+06 -1.0 3.79e-02 - 1.00e+00 1.00e+00h 1 23 2.0000019e+01 1.02e+00 2.16e+06 -1.0 3.75e-02 - 1.00e+00 1.00e+00h 1 24 2.0000018e+01 9.79e-01 2.11e+06 -1.0 3.71e-02 - 1.00e+00 1.00e+00h 1 25 2.0000017e+01 9.25e-01 4.52e+06 -1.7 5.35e-02 - 1.00e+00 1.00e+00h 1 26 2.0000016e+01 8.72e-01 4.41e+06 -1.7 5.30e-02 - 1.00e+00 1.00e+00h 1 27 2.0000016e+01 8.20e-01 4.25e+06 -1.7 5.18e-02 - 1.00e+00 1.00e+00h 1 28 2.0000015e+01 7.70e-01 4.09e+06 -1.7 5.05e-02 - 1.00e+00 1.00e+00h 1 29 2.0000014e+01 7.21e-01 3.94e+06 -1.7 4.91e-02 - 1.00e+00 1.00e+00h 1 iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls 30 2.0000013e+01 6.73e-01 3.78e+06 -1.7 4.77e-02 - 1.00e+00 1.00e+00h 1 31 2.0000013e+01 6.27e-01 3.63e+06 -1.7 4.62e-02 - 1.00e+00 1.00e+00h 1 32 2.0000012e+01 5.98e-01 3.47e+06 -1.7 4.47e-02 - 1.00e+00 1.00e+00h 1 33 2.0000011e+01 5.96e-01 3.32e+06 -1.7 4.32e-02 - 1.00e+00 1.00e+00h 1 34 2.0000011e+01 5.93e-01 3.17e+06 -1.7 4.16e-02 - 1.00e+00 1.00e+00h 1 35 2.0000010e+01 5.91e-01 3.02e+06 -1.7 4.00e-02 - 1.00e+00 1.00e+00h 1 36 2.0000009e+01 5.88e-01 2.88e+06 -1.7 3.84e-02 - 1.00e+00 1.00e+00h 1 37 2.0000009e+01 5.86e-01 2.74e+06 -1.7 3.68e-02 - 1.00e+00 1.00e+00h 1 38 2.0000008e+01 5.83e-01 2.60e+06 -1.7 3.51e-02 - 1.00e+00 1.00e+00h 1 39 2.0000008e+01 5.81e-01 2.47e+06 -1.7 3.35e-02 - 1.00e+00 1.00e+00h 1 iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls 40 2.0000007e+01 5.79e-01 2.34e+06 -1.7 3.20e-02 - 1.00e+00 1.00e+00h 1 41 2.0000007e+01 5.77e-01 2.22e+06 -1.7 3.04e-02 - 1.00e+00 1.00e+00h 1 42 2.0000006e+01 5.74e-01 2.10e+06 -1.7 2.89e-02 - 1.00e+00 1.00e+00h 1 43 2.0000006e+01 5.72e-01 1.98e+06 -1.7 2.74e-02 - 1.00e+00 1.00e+00h 1 44 2.0000005e+01 5.70e-01 1.87e+06 -1.7 2.60e-02 - 1.00e+00 1.00e+00h 1 45 2.0000005e+01 5.68e-01 1.77e+06 -1.7 2.46e-02 - 1.00e+00 1.00e+00h 1 46 2.0000005e+01 5.66e-01 1.67e+06 -1.7 2.33e-02 - 1.00e+00 1.00e+00h 1 47 2.0000004e+01 5.64e-01 1.58e+06 -1.7 2.20e-02 - 1.00e+00 1.00e+00h 1 48 2.0000004e+01 5.62e-01 1.49e+06 -1.7 2.08e-02 - 1.00e+00 1.00e+00h 1 49 2.0000004e+01 5.60e-01 1.40e+06 -1.7 1.96e-02 - 1.00e+00 1.00e+00h 1 iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls 50 2.0000003e+01 5.58e-01 1.32e+06 -1.7 1.85e-02 - 1.00e+00 1.00e+00h 1 51 2.0000003e+01 5.57e-01 1.24e+06 -1.7 1.74e-02 - 1.00e+00 1.00e+00h 1 52 2.0000003e+01 5.55e-01 1.17e+06 -1.7 1.64e-02 - 1.00e+00 1.00e+00h 1 53 2.0000003e+01 5.53e-01 1.10e+06 -1.7 1.55e-02 - 1.00e+00 1.00e+00h 1 54 2.0000002e+01 5.51e-01 1.08e+06 -1.7 1.46e-02 - 1.00e+00 1.00e+00h 1 55 2.0000002e+01 5.50e-01 1.05e+06 -1.7 1.37e-02 - 1.00e+00 1.00e+00h 1 56 2.0000002e+01 5.48e-01 1.03e+06 -1.7 1.29e-02 - 1.00e+00 1.00e+00h 1 57 2.0000002e+01 5.47e-01 1.00e+06 -1.7 1.22e-02 - 1.00e+00 1.00e+00h 1 58 2.0000002e+01 5.45e-01 9.72e+05 -1.7 1.15e-02 - 1.00e+00 1.00e+00h 1 59 2.0000002e+01 5.43e-01 9.45e+05 -1.7 1.08e-02 - 1.00e+00 1.00e+00h 1 iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls 60 2.0000001e+01 5.42e-01 9.17e+05 -1.7 1.02e-02 - 1.00e+00 1.00e+00h 1 61 2.0000001e+01 5.40e-01 8.90e+05 -1.7 9.59e-03 - 1.00e+00 1.00e+00h 1 62 2.0000001e+01 5.39e-01 8.63e+05 -1.7 9.04e-03 - 1.00e+00 1.00e+00h 1 63 2.0000001e+01 5.37e-01 8.37e+05 -1.7 8.52e-03 - 1.00e+00 1.00e+00h 1 64 2.0000001e+01 5.36e-01 8.12e+05 -1.7 8.04e-03 - 1.00e+00 1.00e+00h 1 65 2.0000001e+01 5.35e-01 7.88e+05 -1.7 7.58e-03 - 1.00e+00 1.00e+00h 1 66 2.0000001e+01 5.33e-01 7.64e+05 -1.7 7.15e-03 - 1.00e+00 1.00e+00h 1 67 2.0000001e+01 5.32e-01 7.41e+05 -1.7 6.75e-03 - 1.00e+00 1.00e+00h 1 68 2.0000001e+01 5.30e-01 7.19e+05 -1.7 6.37e-03 - 1.00e+00 1.00e+00h 1 69 2.0000001e+01 5.29e-01 6.98e+05 -1.7 6.02e-03 - 1.00e+00 1.00e+00h 1 iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls 70 2.0000001e+01 5.26e-01 2.63e+06 -1.7 5.69e-03 - 1.00e+00 1.00e+00H 1 71 2.0000000e+01 5.25e-01 6.61e+05 -1.7 5.34e-03 - 1.00e+00 1.00e+00h 1 72 2.0000000e+01 5.23e-01 6.20e+05 -1.7 4.80e-03 - 1.00e+00 1.00e+00h 1 73 2.0000000e+01 5.22e-01 6.01e+05 -1.7 4.54e-03 - 1.00e+00 1.00e+00h 1 74 2.0000000e+01 5.21e-01 5.84e+05 -1.7 4.30e-03 - 1.00e+00 1.00e+00h 1 75 2.0000000e+01 5.19e-01 5.67e+05 -1.7 4.07e-03 - 1.00e+00 1.00e+00h 1 76 2.0000000e+01 5.18e-01 5.51e+05 -1.7 3.86e-03 - 1.00e+00 1.00e+00h 1 77 2.0000000e+01 5.17e-01 5.35e+05 -1.7 3.65e-03 - 1.00e+00 1.00e+00h 1 78 2.0000000e+01 5.14e-01 2.02e+06 -1.7 3.47e-03 - 1.00e+00 1.00e+00H 1 79 2.0000000e+01 5.11e-01 1.94e+06 -1.7 3.30e-03 - 1.00e+00 1.00e+00H 1 iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls 80 2.0000000e+01 5.08e-01 1.83e+06 -1.7 2.99e-03 - 1.00e+00 1.00e+00H 1 81 1.9999997e+01 5.05e-01 1.72e+06 -1.7 2.88e-03 - 1.00e+00 1.00e+00H 1 82 2.0000000e+01 5.04e-01 4.53e+05 -1.7 5.72e-03 - 1.00e+00 1.00e+00s 22 83 2.0000000e+01 5.02e-01 4.44e+05 -1.7 2.74e-03 - 1.00e+00 1.00e+00s 22 84r 2.0000000e+01 5.02e-01 9.99e+02 -0.3 0.00e+00 - 0.00e+00 0.00e+00R 1 85r 2.0000002e+01 6.31e-01 9.93e+02 -0.3 6.70e+01 - 7.42e-03 6.15e-03f 1 86r 2.0000003e+01 6.21e-01 9.89e+02 -0.3 5.36e+00 2.0 1.00e+00 5.36e-03f 1 87r 2.0000005e+01 5.77e-01 8.85e+02 -0.3 8.79e-01 2.4 1.00e+00 1.32e-01f 1 88r 1.9999993e+01 5.81e-01 6.07e+02 -0.3 3.06e-01 2.9 1.00e+00 4.17e-01f 1 89r 1.9999795e+01 5.90e-01 3.23e+02 -0.3 1.70e-01 3.3 1.00e+00 1.00e+00f 1 iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls 90r 1.9999435e+01 5.94e-01 5.80e+01 -0.3 1.03e-02 3.7 1.00e+00 1.00e+00f 1 91r 1.9999237e+01 5.96e-01 5.36e+01 -0.3 3.97e-03 4.1 1.00e+00 1.00e+00f 1 92r 1.9996635e+01 6.00e-01 6.89e+01 -0.3 1.52e-02 3.7 1.00e+00 1.00e+00f 1 93r 1.9995354e+01 6.01e-01 9.57e+01 -0.3 8.03e-03 4.1 1.00e+00 1.00e+00f 1 94r 1.9969366e+01 4.51e-01 6.56e+02 -0.3 1.79e-01 3.6 1.00e+00 1.00e+00f 1 95 1.9957122e+01 1.96e-01 6.25e+07 -1.7 1.68e-01 8.3 1.00e+00 1.00e+00h 1 96 1.9942302e+01 1.25e-01 3.63e+07 -1.7 6.37e-02 8.7 1.00e+00 1.00e+00h 1 97 1.9932539e+01 1.42e-01 2.92e+07 -1.7 2.12e-02 9.2 1.00e+00 1.00e+00h 1 98 1.9893096e+01 1.50e-01 3.36e+07 -1.7 6.96e-02 8.7 1.00e+00 1.00e+00h 1 99 1.9875313e+01 1.39e-01 4.01e+07 -1.7 3.10e-02 9.1 1.00e+00 1.00e+00h 1 iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls 100 1.9809988e+01 8.95e-02 4.42e+07 -1.7 9.54e-02 8.6 1.00e+00 1.00e+00h 1 101 1.9780147e+01 6.48e-02 5.55e+07 -1.7 4.84e-02 9.1 1.00e+00 1.00e+00h 1 102 1.9632464e+01 5.57e-01 1.17e+08 -1.7 3.37e-01 8.6 1.00e+00 1.00e+00h 1 103 1.9561605e+01 3.18e-01 9.34e+07 -1.7 1.15e-01 9.0 1.00e+00 1.00e+00h 1 104 1.9383371e+01 2.49e-01 4.98e+07 -1.7 1.47e-01 8.5 1.00e+00 1.00e+00h 1 105 1.9212168e+01 2.85e-01 1.78e+07 -1.7 1.28e-01 8.1 1.00e+00 1.00e+00h 1 106 1.9276108e+01 2.13e-01 2.19e+07 -1.7 7.55e-02 8.5 1.00e+00 1.00e+00h 1 107 1.9574712e+01 1.48e-01 2.74e+07 -1.7 2.90e-01 8.0 1.00e+00 1.00e+00h 1 108 1.9677940e+01 5.86e-02 2.13e+07 -1.7 8.35e-02 8.4 1.00e+00 1.00e+00h 1 109 1.9673235e+01 2.30e-02 5.12e+06 -1.7 1.69e-01 8.0 1.00e+00 1.00e+00H 1 iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls 110 1.9672936e+01 4.46e-02 1.07e+06 -1.7 5.17e-02 7.5 1.00e+00 1.00e+00h 1 111 1.9665062e+01 1.86e-02 4.37e+05 -1.7 4.55e-02 7.0 1.00e+00 1.00e+00h 1 112 1.9639714e+01 7.62e-03 9.64e+04 -1.7 3.78e-02 6.5 1.00e+00 1.00e+00h 1 113 1.9643614e+01 1.02e-03 5.07e+03 -1.7 6.16e-03 6.0 1.00e+00 1.00e+00h 1 114 1.9644944e+01 5.14e-05 6.85e+02 -1.7 1.82e-03 5.6 1.00e+00 1.00e+00h 1 115 1.9644828e+01 6.86e-06 9.14e+00 -1.7 7.38e-05 5.1 1.00e+00 1.00e+00h 1 116 1.9644635e+01 8.62e-08 1.13e+00 -1.7 2.76e-05 4.6 1.00e+00 1.00e+00h 1 117 1.9644081e+01 4.76e-09 1.10e+00 -1.7 8.08e-05 4.1 1.00e+00 1.00e+00h 1 118 1.9642421e+01 6.37e-08 1.10e+00 -1.7 2.42e-04 3.7 1.00e+00 1.00e+00f 1 119 1.9637439e+01 5.73e-07 1.10e+00 -1.7 7.27e-04 3.2 1.00e+00 1.00e+00f 1 iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls 120 1.9622487e+01 5.16e-06 1.10e+00 -1.7 2.18e-03 2.7 1.00e+00 1.00e+00f 1 121 1.9577559e+01 4.64e-05 1.11e+00 -1.7 6.55e-03 2.2 1.00e+00 1.00e+00f 1 122 1.9442160e+01 4.19e-04 1.11e+00 -1.7 1.97e-02 1.8 1.00e+00 1.00e+00f 1 123 1.9030632e+01 3.78e-03 1.11e+00 -1.7 5.97e-02 1.3 1.00e+00 1.00e+00f 1 124 1.7754603e+01 3.33e-02 1.08e+00 -1.7 1.80e-01 0.8 1.00e+00 1.00e+00f 1 125 1.3741171e+01 2.93e-01 1.11e+00 -1.7 5.35e-01 0.3 1.00e+00 1.00e+00f 1 126 2.3733321e+00 3.49e+00 1.22e+00 -1.7 1.69e+00 -0.2 1.00e+00 1.00e+00f 1 127 2.2985871e+00 2.25e+00 1.70e+00 -1.7 8.91e-01 0.3 1.00e+00 1.00e+00h 1 128 2.6257083e+00 1.08e+00 2.12e+00 -1.7 6.37e-01 0.7 1.00e+00 1.00e+00h 1 129 4.9342499e+00 5.74e-01 9.46e-01 -1.7 6.42e-01 - 1.00e+00 1.00e+00h 1 iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls 130 6.9097741e+00 6.91e-02 2.42e-01 -1.7 5.72e-01 - 1.00e+00 1.00e+00h 1 131 7.3451482e+00 3.29e-03 4.20e-02 -1.7 1.15e-01 - 1.00e+00 1.00e+00h 1 132 7.3729165e+00 5.26e-05 3.16e-04 -2.5 8.96e-03 - 1.00e+00 1.00e+00h 1 133 7.3730955e+00 7.89e-09 3.92e-08 -5.7 8.63e-05 - 1.00e+00 1.00e+00h 1 134 7.3730955e+00 4.00e-15 1.78e-15 -8.6 8.04e-09 - 1.00e+00 1.00e+00h 1 Number of Iterations....: 134 (scaled) (unscaled) Objective...............: 7.3730955000280947e+00 7.3730955000280947e+00 Dual infeasibility......: 1.7763568394002505e-15 1.7763568394002505e-15 Constraint violation....: 3.9968028886505635e-15 3.9968028886505635e-15 Complementarity.........: 0.0000000000000000e+00 0.0000000000000000e+00 Overall NLP error.......: 3.9968028886505635e-15 3.9968028886505635e-15 Number of objective function evaluations = 166 Number of objective gradient evaluations = 126 Number of equality constraint evaluations = 166 Number of inequality constraint evaluations = 0 Number of equality constraint Jacobian evaluations = 136 Number of inequality constraint Jacobian evaluations = 0 Number of Lagrangian Hessian evaluations = 134 Total CPU secs in IPOPT (w/o function evaluations) = 0.410 Total CPU secs in NLP function evaluations = 0.116 EXIT: Optimal Solution Found.
solution_summary(model)
* Solver : Ipopt * Status Termination status : LOCALLY_SOLVED Primal status : FEASIBLE_POINT Dual status : FEASIBLE_POINT Message from the solver: "Solve_Succeeded" * Candidate solution Objective value : 7.373095500028095 * Work counters Solve time (sec) : 0.57603
Pojďme se výslednou křivku nakreslit!
using PyPlot
Souřadnice bodů:
xs = value.(x)
ys = value.(y);
A graf:
fig, ax = plt.subplots()
ax.axis("equal")
ax.grid()
scatter([0, D], [a, a], 50, color="red")
plot(xs, ys)
1-element Vector{PyCall.PyObject}: PyObject <matplotlib.lines.Line2D object at 0x7f64c0419670>
Tato křivka se dá popsat analyticky pomocí kosinu hyperbolického a má dokonce jméno, jde o catenarii, česky "řetězovku". Sice připomíná parabolu, ale není to parabola!
Dále k této lekci přikládáme Pluto notebook, kde si s touto úlohou můžete interaktivněji pohrát.
3.3 Sudoku
Populární hlavolam Sudoku asi není potřeba představovat.
Zde si ukážeme, jak problém nalezení řešení Sudoku hlavolamu lze zformulovat jakožto lineární program s binárními proměnnými. Objektivní funkci bychom mohli formálně zadat třeba konstantní, není potřeba. Budeme hledat pouze nějaké přípustné řešení.
Použijeme GLPK řešič, který umí řešit úlohy lineární programování i s celočíselnými proměnnými (binární proměnné do této škatulky také spadnou).
Nezapomeňte nainstalovat Julia balíček poskytující tento řešič, tj. ] add GLPK
.
using JuMP, GLPK
Klíčem k úspěchu je následující pohled na sudoku problém: místo čísel až v matici si představme jedničky a nuly v krychli . Krychle je vyplněna nulami, jen a pouze je-li v původní matici hodnota , pak tomu odpovídá hodnota .
Nyní se stačí zamyslet jak popsat podmínky na správné Sudoku řešení (cifry se nesmí opakovat v řádcích, sloupcích a v devíti podblocích. Zároveň musíme ctít zadání.
function sudoku_solver(puzzle::Matrix{Int64})
if size(puzzle, 1) != 9 || size(puzzle, 2) != 9
error("Sudoku puzzle has to be a 9x9 matrix!")
end
# let's use the GLPK solver
model = Model(GLPK.Optimizer)
# now we define our variables, 9^3 of them, do not be afraid!
@variable(model, x[i=1:9, j=1:9, k=1:9], Bin)
for i in 1:9, j in 1:9
# only one digit at each site
@constraint(model, sum(x[i, j, k] for k=1:9) == 1)
# fix prescribed variables
if puzzle[i, j] != 0
fix(x[i, j, puzzle[i, j]], 1, force=true)
end
end
for i=1:9, k=1:9
@constraint(model, sum(x[i, j, k] for j=1:9) == 1)
@constraint(model, sum(x[j, i, k] for j=1:9) == 1)
end
for i=0:2, j=0:2, k=1:9
@constraint(model, sum(x[3i+l, 3j+m, k] for l=1:3, m=1:3) == 1)
end
# solve it!
optimize!(model)
@info solution_summary(model)
# extract a solution a return it in the matrix form
solution = value.(x)
return Int64[ findfirst(isone, solution[i, j, :]) for i=1:9, j=1:9 ]
end
sudoku_solver (generic function with 1 method)
Zkusme vyřešit nějaké Sudoku a nebuďme přízemní. Následující zadání je Evil obtížnost z webu sudoku.com:
Sudoku evil is a 9×9 grid sudoku puzzle with the highest possible level of difficulty. It’s an advanced level only for experienced sudoku solvers. Simple logic and basic knowledge of sudoku rules won’t cut it – you as a player should know advanced sudoku solving techniques and understand how to apply them in practice.
🤘 = [
4 0 0 0 7 8 0 0 3
0 0 0 2 0 0 0 0 0
0 0 9 0 0 0 0 1 0
0 1 0 0 6 2 0 3 0
0 0 0 4 0 0 0 0 2
6 0 0 5 0 0 0 0 0
0 4 0 0 0 0 7 0 0
7 0 0 0 3 6 0 0 8
0 0 0 0 5 0 0 0 0
];
sudoku_solver(🤘)
┌ Info: * Solver : GLPK │ │ * Status │ Termination status : OPTIMAL │ Primal status : FEASIBLE_POINT │ Dual status : NO_SOLUTION │ Message from the solver: │ "Solution is optimal" │ │ * Candidate solution │ Objective value : 0.0 │ Objective bound : 0.0 │ │ * Work counters │ Solve time (sec) : 0.03582 └ @ Main In[24]:33
9×9 Matrix{Int64}: 4 5 1 6 7 8 9 2 3 3 7 6 2 9 1 4 8 5 2 8 9 3 4 5 6 1 7 5 1 7 9 6 2 8 3 4 9 3 8 4 1 7 5 6 2 6 2 4 5 8 3 1 7 9 1 4 3 8 2 9 7 5 6 7 9 5 1 3 6 2 4 8 8 6 2 7 5 4 3 9 1
3.4 Exponenciální fit
V poslední době jste jistě zaslechli o exponenciálním růstu, resp. funkci, a odhadu parametrů tohoto růstu z dat. Jak v časové řadě (posloupnost číselných dat) takovéto chování odhalit a jak odhadnout jeho parametry?
K tomu lze využít metodu nejmenších čtverců (není to jediný způsob). Předpokládejme, že máme naměřené hodnoty nějaké veličiny v časových okamžicích (nezávisle proměnná nutně nemusí být čas ). Dále máme podezření, že hypotetická závislost se chová exponenciálně, obecně třeba ve tvaru , kde je trojice neznámých reálných parametrů.
Naším cílem je nalézt hodnoty těchto neznámých parametrů tak, aby výsledná funkční závislost co nejlépe vystihovala naměřená data. Ve smyslu metody nejmenších čtverců tím máme na mysli, aby funkce nabývala své minimální hodnoty (tj. aby součet kvadrátů odchylek byl nejmenší možný, odtud název metody).
Očividně se jedná o nelineární problém. Použijeme Ipopt řešič (] add Ipopt
).
using JuMP, Ipopt, PyPlot
function expfit(ts, ys; max_iter=100, optimizer=Ipopt.Optimizer)
model = Model(optimizer)
set_optimizer_attribute(model, "max_iter", max_iter)
if max_iter > 100
set_silent(model)
end
# variables
@variable(model, a[i = 1:3])
# non-linear objective function
@NLobjective(model, Min, sum((a[1] * exp(a[2] * ts[j]) + a[3] - ys[j])^2 for j=1:length(ts)))
# solve it!
optimize!(model)
@info solution_summary(model)
return objective_value(model), value.(a)
end
expfit (generic function with 1 method)
Příklad č. 1
Testovací data, exponenciální růst plus šum.
ts = LinRange(-1, 4, 50)
ys = [ 0.323 * exp(1.56 * t) - 2.453 + 4*(rand() - 0.5) for t in ts ];
Grafické znázornění.
fig, ax = plt.subplots()
ax.grid()
ax.set_xlabel(raw"$t$")
ax.set_ylabel(raw"$y$")
scatter(ts, ys, s=2);
_, params = expfit(ts, ys)
****************************************************************************** This program contains Ipopt, a library for large-scale nonlinear optimization. Ipopt is released as open source code under the Eclipse Public License (EPL). For more information visit https://github.com/coin-or/Ipopt ****************************************************************************** This is Ipopt version 3.13.4, running with linear solver mumps. NOTE: Other linear solvers might be more efficient (see Ipopt documentation). Number of nonzeros in equality constraint Jacobian...: 0 Number of nonzeros in inequality constraint Jacobian.: 0 Number of nonzeros in Lagrangian Hessian.............: 6 Total number of variables............................: 3 variables with only lower bounds: 0 variables with lower and upper bounds: 0 variables with only upper bounds: 0 Total number of equality constraints.................: 0 Total number of inequality constraints...............: 0 inequality constraints with only lower bounds: 0 inequality constraints with lower and upper bounds: 0 inequality constraints with only upper bounds: 0 iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls 0 9.5227720e+04 0.00e+00 1.00e+02 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0 1 9.5187171e+04 0.00e+00 1.00e+02 -1.0 1.00e-02 4.0 1.00e+00 1.00e+00f 1 2 9.5063826e+04 0.00e+00 1.01e+02 -1.0 3.04e-02 3.5 1.00e+00 1.00e+00f 1 3 9.4603734e+04 0.00e+00 1.21e+02 -1.0 1.07e-01 3.0 1.00e+00 1.00e+00f 1 4 9.4370065e+04 0.00e+00 1.35e+02 -1.0 4.53e-02 3.5 1.00e+00 1.00e+00f 1 5 8.8923633e+04 0.00e+00 1.02e+03 -1.0 3.92e-01 3.0 1.00e+00 1.00e+00f 1 6 8.7092700e+04 0.00e+00 1.33e+03 -1.0 6.14e-02 4.3 1.00e+00 1.00e+00f 1 7 4.5427446e+04 0.00e+00 7.28e+03 -1.0 9.39e-01 3.8 1.00e+00 5.00e-01f 2 8 2.2442352e+04 0.00e+00 2.39e+04 -1.0 3.05e+00 4.3 1.00e+00 1.25e-01f 4 9 4.7112071e+03 0.00e+00 7.53e+03 -1.0 5.09e-02 4.7 1.00e+00 1.00e+00f 1 iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls 10 1.6139004e+03 0.00e+00 1.56e+03 -1.0 2.74e-02 4.2 1.00e+00 1.00e+00f 1 11 1.3554412e+03 0.00e+00 1.62e+02 -1.0 2.82e-02 3.7 1.00e+00 1.00e+00f 1 12 1.1738488e+03 0.00e+00 1.64e+02 -1.0 6.58e-02 3.3 1.00e+00 1.00e+00f 1 13 9.0920684e+02 0.00e+00 2.01e+03 -1.0 8.97e-01 - 1.00e+00 2.50e-01f 3 14 1.1735797e+02 0.00e+00 4.94e+02 -1.0 3.49e+00 - 1.00e+00 1.00e+00f 1 15 8.8001142e+01 0.00e+00 1.59e+02 -1.0 4.32e-02 2.8 1.00e+00 1.00e+00f 1 16 6.6444342e+01 0.00e+00 2.03e+02 -1.0 7.80e-01 - 1.00e+00 1.00e+00f 1 17 6.1418292e+01 0.00e+00 2.34e+01 -1.0 1.10e-01 - 1.00e+00 1.00e+00f 1 18 6.0720789e+01 0.00e+00 2.25e+01 -1.0 9.79e-02 - 1.00e+00 1.00e+00f 1 19 6.0679659e+01 0.00e+00 1.04e-01 -1.0 8.24e-03 - 1.00e+00 1.00e+00f 1 iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls 20 6.0679609e+01 0.00e+00 1.86e-03 -2.5 8.88e-04 - 1.00e+00 1.00e+00f 1 21 6.0679609e+01 0.00e+00 1.29e-10 -3.8 2.56e-07 - 1.00e+00 1.00e+00f 1 Number of Iterations....: 21 (scaled) (unscaled) Objective...............: 2.9924252062699130e+00 6.0679609324508796e+01 Dual infeasibility......: 1.2924596324375634e-10 2.6208155645690567e-09 Constraint violation....: 0.0000000000000000e+00 0.0000000000000000e+00 Complementarity.........: 0.0000000000000000e+00 0.0000000000000000e+00 Overall NLP error.......: 1.2924596324375634e-10 2.6208155645690567e-09 Number of objective function evaluations = 40 Number of objective gradient evaluations = 22 Number of equality constraint evaluations = 0 Number of inequality constraint evaluations = 0 Number of equality constraint Jacobian evaluations = 0 Number of inequality constraint Jacobian evaluations = 0 Number of Lagrangian Hessian evaluations = 21 Total CPU secs in IPOPT (w/o function evaluations) = 0.875 Total CPU secs in NLP function evaluations = 1.408 EXIT: Optimal Solution Found.
┌ Info: * Solver : Ipopt │ │ * Status │ Termination status : LOCALLY_SOLVED │ Primal status : FEASIBLE_POINT │ Dual status : FEASIBLE_POINT │ Message from the solver: │ "Solve_Succeeded" │ │ * Candidate solution │ Objective value : 60.679609324508796 │ │ * Work counters │ Solve time (sec) : 2.60207 └ @ Main In[29]:17
(60.679609324508796, [0.3666528949394605, 1.5257033871545445, -2.4547198128810988])
Prostým porovnáním čísel, která jsme k vygenerování dat použili je vidět, že jsme dostali asi docela rozumné řešení. Znázorněme ho graficky.
fig, ax = plt.subplots()
ax.grid()
ax.set_xlabel(raw"$t$")
ax.set_ylabel(raw"$y$")
scatter(ts, ys, s=3);
plot(ts, [ params[1] * exp(params[2] * t) + params[3] for t in ts ], "r");
Příklad č. 2
Testovací data, exponenciála nemusí nutně jenom růst!
ts = LinRange(0, 7, 100)
ys = [ -2.111 * exp(-0.687 * t) + 8.453 + (rand() - 0.5) / (t^2 + 1) for t in ts ];
Grafické znázornění.
fig, ax = plt.subplots()
ax.grid()
ax.set_xlabel(raw"$t$")
ax.set_ylabel(raw"$y$")
scatter(ts, ys, s=2);
_, params = expfit(ts, ys)
This is Ipopt version 3.13.4, running with linear solver mumps. NOTE: Other linear solvers might be more efficient (see Ipopt documentation). Number of nonzeros in equality constraint Jacobian...: 0 Number of nonzeros in inequality constraint Jacobian.: 0 Number of nonzeros in Lagrangian Hessian.............: 6 Total number of variables............................: 3 variables with only lower bounds: 0 variables with lower and upper bounds: 0 variables with only upper bounds: 0 Total number of equality constraints.................: 0 Total number of inequality constraints...............: 0 inequality constraints with only lower bounds: 0 inequality constraints with lower and upper bounds: 0 inequality constraints with only upper bounds: 0 iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls 0 6.4335259e+03 0.00e+00 1.00e+02 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0 1 6.4015984e+03 0.00e+00 9.99e+01 -1.0 9.99e-03 4.0 1.00e+00 1.00e+00f 1 2 6.3052332e+03 0.00e+00 1.01e+02 -1.0 3.02e-02 3.5 1.00e+00 1.00e+00f 1 3 5.9541964e+03 0.00e+00 1.18e+02 -1.0 1.05e-01 3.0 1.00e+00 1.00e+00f 1 4 5.7740252e+03 0.00e+00 1.32e+02 -1.0 4.43e-02 3.5 1.00e+00 1.00e+00f 1 5 3.2789101e+03 0.00e+00 5.60e+02 -1.0 6.22e-01 3.0 1.00e+00 5.00e-01f 2 6 3.0576021e+03 0.00e+00 1.35e+03 -1.0 5.35e-01 2.5 1.00e+00 2.50e-01f 3 7 1.9904464e+03 0.00e+00 1.34e+02 -1.0 3.89e-01 2.9 1.00e+00 1.00e+00f 1 8 1.0632243e+02 0.00e+00 1.30e+02 -1.0 6.92e+00 - 1.00e+00 1.00e+00f 1 9 3.5323668e+01 0.00e+00 4.58e+01 -1.0 4.12e-02 2.5 1.00e+00 1.00e+00f 1 iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls 10 1.6791890e+01 0.00e+00 1.34e+01 -1.0 3.11e-02 2.0 1.00e+00 1.00e+00f 1 11 1.4013027e+01 0.00e+00 2.62e+00 -1.0 2.13e-02 1.5 1.00e+00 1.00e+00f 1 12 1.3770453e+01 0.00e+00 1.64e-01 -1.0 1.92e-02 1.0 1.00e+00 1.00e+00f 1 13 1.3594999e+01 0.00e+00 6.94e+00 -2.5 5.82e-01 - 1.00e+00 5.00e-01f 2 14 1.2517729e+01 0.00e+00 2.95e+00 -2.5 2.00e-01 - 1.00e+00 1.00e+00f 1 15 1.2229710e+01 0.00e+00 1.58e-01 -2.5 4.49e-02 0.6 1.00e+00 1.00e+00f 1 16 1.1969189e+01 0.00e+00 6.37e+00 -2.5 6.95e-01 - 1.00e+00 5.00e-01f 2 17 1.1350231e+01 0.00e+00 6.22e-01 -2.5 2.20e-01 - 1.00e+00 1.00e+00f 1 18 1.1098839e+01 0.00e+00 5.61e+00 -2.5 1.59e+00 - 1.00e+00 2.50e-01f 3 19 1.0722710e+01 0.00e+00 5.92e-01 -2.5 2.61e-01 - 1.00e+00 1.00e+00f 1 iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls 20 1.0540648e+01 0.00e+00 5.36e+00 -2.5 8.95e-01 - 1.00e+00 5.00e-01f 2 21 1.0282819e+01 0.00e+00 1.20e+00 -2.5 3.60e-01 - 1.00e+00 1.00e+00f 1 22 1.0145434e+01 0.00e+00 4.29e+00 -2.5 9.16e-01 - 1.00e+00 5.00e-01f 2 23 9.9696052e+00 0.00e+00 2.90e+00 -2.5 5.16e-01 - 1.00e+00 1.00e+00f 1 24 9.8449345e+00 0.00e+00 6.62e+00 -2.5 7.42e-01 - 1.00e+00 1.00e+00f 1 25 9.7021769e+00 0.00e+00 1.47e+00 -2.5 5.16e-01 - 1.00e+00 1.00e+00f 1 26 9.6269076e+00 0.00e+00 5.55e+00 -2.5 1.39e+00 - 1.00e+00 5.00e-01f 2 27 9.5270982e+00 0.00e+00 3.02e+00 -2.5 7.00e-01 - 1.00e+00 1.00e+00f 1 28 9.4705872e+00 0.00e+00 1.00e+01 -2.5 1.15e+00 - 1.00e+00 1.00e+00f 1 29 9.3748495e+00 0.00e+00 7.94e-01 -2.5 6.14e-01 - 1.00e+00 1.00e+00f 1 iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls 30 9.3553504e+00 0.00e+00 1.34e+01 -2.5 2.81e+00 - 1.00e+00 5.00e-01f 2 31 9.2642589e+00 0.00e+00 2.06e-01 -2.5 6.31e-01 - 1.00e+00 1.00e+00f 1 32 9.2304760e+00 0.00e+00 7.23e+00 -2.5 4.70e+00 - 1.00e+00 2.50e-01f 3 33 9.1850111e+00 0.00e+00 3.66e+00 -2.5 1.06e+00 - 1.00e+00 1.00e+00f 1 34 9.1639452e+00 0.00e+00 1.39e+01 -2.5 1.87e+00 - 1.00e+00 1.00e+00f 1 35 9.1138112e+00 0.00e+00 7.48e-01 -2.5 8.45e-01 - 1.00e+00 1.00e+00f 1 36 9.0935234e+00 0.00e+00 6.18e+00 -2.5 5.16e+00 - 1.00e+00 2.50e-01f 3 37 9.0672335e+00 0.00e+00 7.38e+00 -2.5 1.66e+00 - 1.00e+00 1.00e+00f 1 38 9.0426079e+00 0.00e+00 5.80e+00 -2.5 1.60e+00 - 1.00e+00 1.00e+00f 1 39 9.0231881e+00 0.00e+00 9.95e+00 -2.5 2.07e+00 - 1.00e+00 1.00e+00f 1 iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls 40 9.0019896e+00 0.00e+00 3.81e+00 -2.5 1.51e+00 - 1.00e+00 1.00e+00f 1 41 8.9901446e+00 0.00e+00 7.51e+00 -2.5 3.29e+00 - 1.00e+00 5.00e-01f 2 42 8.9742011e+00 0.00e+00 7.84e+00 -2.5 2.12e+00 - 1.00e+00 1.00e+00f 1 43 8.9597379e+00 0.00e+00 7.86e+00 -2.5 2.21e+00 - 1.00e+00 1.00e+00f 1 44 8.9467893e+00 0.00e+00 8.50e+00 -2.5 2.38e+00 - 1.00e+00 1.00e+00f 1 45 8.9348260e+00 0.00e+00 7.93e+00 -2.5 2.41e+00 - 1.00e+00 1.00e+00f 1 46 8.9243786e+00 0.00e+00 9.78e+00 -2.5 2.72e+00 - 1.00e+00 1.00e+00f 1 47 8.9142045e+00 0.00e+00 7.06e+00 -2.5 2.48e+00 - 1.00e+00 1.00e+00f 1 48 8.9065820e+00 0.00e+00 1.41e+01 -2.5 3.42e+00 - 1.00e+00 1.00e+00f 1 49 8.8970995e+00 0.00e+00 3.87e+00 -2.5 2.14e+00 - 1.00e+00 1.00e+00f 1 iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls 50 8.8919516e+00 0.00e+00 1.24e+01 -2.5 6.15e+00 - 1.00e+00 5.00e-01f 2 51 8.8843562e+00 0.00e+00 5.82e+00 -2.5 2.64e+00 - 1.00e+00 1.00e+00f 1 52 8.8824205e+00 0.00e+00 2.55e+01 -2.5 5.09e+00 - 1.00e+00 1.00e+00f 1 53 8.8716934e+00 0.00e+00 1.04e+00 -2.5 1.76e+00 - 1.00e+00 1.00e+00f 1 54 8.8682122e+00 0.00e+00 1.64e+01 -2.5 1.67e+01 - 1.00e+00 2.50e-01f 3 55 8.8622288e+00 0.00e+00 4.23e+00 -2.5 2.66e+00 - 1.00e+00 1.00e+00f 1 56 8.8590651e+00 0.00e+00 1.51e+01 -2.5 8.21e+00 - 1.00e+00 5.00e-01f 2 57 8.8541282e+00 0.00e+00 5.66e+00 -2.5 3.13e+00 - 1.00e+00 1.00e+00f 1 58 8.8513011e+00 0.00e+00 1.19e+01 -2.5 7.30e+00 - 1.00e+00 5.00e-01f 2 59 8.8473437e+00 0.00e+00 1.03e+01 -2.5 4.19e+00 - 1.00e+00 1.00e+00f 1 iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls 60 8.8439356e+00 0.00e+00 1.43e+01 -2.5 4.96e+00 - 1.00e+00 1.00e+00f 1 61 8.8403537e+00 0.00e+00 7.98e+00 -2.5 3.99e+00 - 1.00e+00 1.00e+00f 1 62 8.8385753e+00 0.00e+00 2.56e+01 -2.5 6.81e+00 - 1.00e+00 1.00e+00f 1 63 8.8342672e+00 0.00e+00 2.45e+00 -2.5 2.79e+00 - 1.00e+00 1.00e+00f 1 64 8.8324284e+00 0.00e+00 1.22e+01 -2.5 1.81e+01 - 1.00e+00 2.50e-01f 3 65 8.8298738e+00 0.00e+00 1.34e+01 -2.5 5.50e+00 - 1.00e+00 1.00e+00f 1 66 8.8274105e+00 0.00e+00 1.18e+01 -2.5 5.34e+00 - 1.00e+00 1.00e+00f 1 67 8.8252595e+00 0.00e+00 1.60e+01 -2.5 6.25e+00 - 1.00e+00 1.00e+00f 1 68 8.8230015e+00 0.00e+00 9.33e+00 -2.5 5.09e+00 - 1.00e+00 1.00e+00f 1 69 8.8217571e+00 0.00e+00 2.73e+01 -2.5 8.38e+00 - 1.00e+00 1.00e+00f 1 iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls 70 8.8191175e+00 0.00e+00 3.22e+00 -2.5 3.60e+00 - 1.00e+00 1.00e+00f 1 71 8.8179646e+00 0.00e+00 1.24e+01 -2.5 2.10e+01 - 1.00e+00 2.50e-01f 3 72 8.8163997e+00 0.00e+00 1.83e+01 -2.5 7.46e+00 - 1.00e+00 1.00e+00f 1 73 8.8146783e+00 0.00e+00 8.94e+00 -2.5 5.60e+00 - 1.00e+00 1.00e+00f 1 74 8.8142495e+00 0.00e+00 3.68e+01 -2.5 1.08e+01 - 1.00e+00 1.00e+00f 1 75 8.8116978e+00 0.00e+00 1.97e+00 -2.5 3.49e+00 - 1.00e+00 1.00e+00f 1 76 8.8108983e+00 0.00e+00 2.50e+01 -2.5 3.57e+01 - 1.00e+00 2.50e-01f 3 77 8.8093552e+00 0.00e+00 5.60e+00 -2.5 5.07e+00 - 1.00e+00 1.00e+00f 1 78 8.8086497e+00 0.00e+00 2.64e+01 -2.5 1.87e+01 - 1.00e+00 5.00e-01f 2 79 8.8072332e+00 0.00e+00 5.46e+00 -2.5 5.22e+00 - 1.00e+00 1.00e+00f 1 iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls 80 8.8066524e+00 0.00e+00 2.93e+01 -2.5 2.06e+01 - 1.00e+00 5.00e-01f 2 81 8.8052956e+00 0.00e+00 4.72e+00 -2.5 5.15e+00 - 1.00e+00 1.00e+00f 1 82 8.8050445e+00 0.00e+00 3.85e+01 -2.5 2.48e+01 - 1.00e+00 5.00e-01f 2 83 8.8034728e+00 0.00e+00 2.68e+00 -2.5 4.44e+00 - 1.00e+00 1.00e+00f 1 84 8.8028870e+00 0.00e+00 2.36e+01 -2.5 4.00e+01 - 1.00e+00 2.50e-01f 3 85 8.8019328e+00 0.00e+00 8.85e+00 -2.5 7.10e+00 - 1.00e+00 1.00e+00f 1 86 8.8013765e+00 0.00e+00 1.89e+01 -2.5 1.72e+01 - 1.00e+00 5.00e-01f 2 87 8.8005720e+00 0.00e+00 1.48e+01 -2.5 9.15e+00 - 1.00e+00 1.00e+00f 1 88 8.7999148e+00 0.00e+00 2.49e+01 -2.5 1.18e+01 - 1.00e+00 1.00e+00f 1 89 8.7991219e+00 0.00e+00 9.21e+00 -2.5 7.75e+00 - 1.00e+00 1.00e+00f 1 iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls 90 8.7986598e+00 0.00e+00 2.02e+01 -2.5 1.92e+01 - 1.00e+00 5.00e-01f 2 91 8.7979868e+00 0.00e+00 1.50e+01 -2.5 9.90e+00 - 1.00e+00 1.00e+00f 1 92 8.7974587e+00 0.00e+00 2.77e+01 -2.5 1.33e+01 - 1.00e+00 1.00e+00f 1 93 8.7967777e+00 0.00e+00 8.53e+00 -2.5 8.05e+00 - 1.00e+00 1.00e+00f 1 94 8.7963967e+00 0.00e+00 2.47e+01 -2.5 2.34e+01 - 1.00e+00 5.00e-01f 2 95 8.7958042e+00 0.00e+00 1.16e+01 -2.5 9.45e+00 - 1.00e+00 1.00e+00f 1 96 8.7957438e+00 0.00e+00 5.19e+01 -2.5 1.91e+01 - 1.00e+00 1.00e+00f 1 97 8.7947533e+00 0.00e+00 2.34e+00 -2.5 5.44e+00 - 1.00e+00 1.00e+00f 1 98 8.7945912e+00 0.00e+00 4.34e+01 -2.5 7.03e+01 - 1.00e+00 2.50e-01f 3 99 8.7938833e+00 0.00e+00 3.98e+00 -2.5 6.64e+00 - 1.00e+00 1.00e+00f 1 iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls 100 8.7935882e+00 0.00e+00 2.41e+01 -2.5 5.21e+01 - 1.00e+00 2.50e-01f 3 Number of Iterations....: 100 (scaled) (unscaled) Objective...............: 5.4950736448569415e-01 8.7935882183300187e+00 Dual infeasibility......: 2.4123397326357445e+01 3.8603890725594329e+02 Constraint violation....: 0.0000000000000000e+00 0.0000000000000000e+00 Complementarity.........: 0.0000000000000000e+00 0.0000000000000000e+00 Overall NLP error.......: 2.4123397326357445e+01 3.8603890725594329e+02 Number of objective function evaluations = 252 Number of objective gradient evaluations = 101 Number of equality constraint evaluations = 0 Number of inequality constraint evaluations = 0 Number of equality constraint Jacobian evaluations = 0 Number of inequality constraint Jacobian evaluations = 0 Number of Lagrangian Hessian evaluations = 100 Total CPU secs in IPOPT (w/o function evaluations) = 0.070 Total CPU secs in NLP function evaluations = 0.007 EXIT: Maximum Number of Iterations Exceeded.
┌ Info: * Solver : Ipopt │ │ * Status │ Termination status : ITERATION_LIMIT │ Primal status : UNKNOWN_RESULT_STATUS │ Dual status : UNKNOWN_RESULT_STATUS │ Message from the solver: │ "Maximum_Iterations_Exceeded" │ │ * Candidate solution │ Objective value : 8.793588218330019 │ │ * Work counters │ Solve time (sec) : 0.08143 └ @ Main In[29]:17
(8.793588218330019, [408.0376224666501, 0.0005687800980638853, -400.85069728568226])
Znázorněme ho graficky.
fig, ax = plt.subplots()
ax.grid()
ax.set_xlabel(raw"$t$")
ax.set_ylabel(raw"$y$")
scatter(ts, ys, s=3);
plot(ts, [ params[1] * exp(params[2] * t) + params[3] for t in ts ], "r");
To nevypadá příliš dobře, dejme řešiči více možnost více iterací.
_, params = expfit(ts, ys, max_iter=10_000)
┌ Info: * Solver : Ipopt │ │ * Status │ Termination status : ITERATION_LIMIT │ Primal status : UNKNOWN_RESULT_STATUS │ Dual status : UNKNOWN_RESULT_STATUS │ Message from the solver: │ "Maximum_Iterations_Exceeded" │ │ * Candidate solution │ Objective value : 8.780490107554508 │ │ * Work counters │ Solve time (sec) : 19.16486 └ @ Main In[29]:17
(8.780490107554508, [-3.115379883761534e6, -7.47716205890407e-8, 3.115387069802581e6])
fig, ax = plt.subplots()
ax.grid()
ax.set_xlabel(raw"$t$")
ax.set_ylabel(raw"$y$")
scatter(ts, ys, s=3);
plot(ts, [ params[1] * exp(params[2] * t) + params[3] for t in ts ], "r");
Chování Ipopt je v tomto případě lehce náhodné!
Řešení některých úloh
"Zemědělské" podmínky.
@variable(model, x[1:3] >= 0)
@constraint(model, x[1] + x[2] + x[3] == S)
@objective(model, Max, (v1 - p1 - o1) * x[1] + (v2 - p2 - o2) * x[2] - (v3 - p3 - o3) * x[3])
Sudoku solver.
function sudoku_solver(puzzle::Matrix{Int64})
if size(puzzle, 1) != 9 || size(puzzle, 2) != 9
error("Sudoku puzzle has to be a 9x9 matrix!")
end
# let's use the GLPK solver
model = Model(GLPK.Optimizer)
# now we define our variables, 9^3 of them, do not be afraid!
@variable(model, x[i=1:9, j=1:9, k=1:9], Bin)
for i in 1:9, j in 1:9
# only one digit at each site
@constraint(model, sum(x[i, j, k] for k=1:9) == 1)
# fix prescribed variables
if puzzle[i, j] != 0
fix(x[i, j, puzzle[i, j]], 1, force=true)
end
end
for i in 1:9, k in 1:9
# constraint rows
@constraint(model, sum(x[i, j, k] for j=1:9) == 1)
# constraint columns
@constraint(model, sum(x[j, i, k] for j=1:9) == 1)
end
# constraint boxes
for is in (1,4,7), js in (1,4,7), k in 1:9
@constraint(model, sum(x[i, j, k] for i=is:(is+2), j=js:(js+2)) == 1)
end
# solve it!
optimize!(model)
# @info solution_summary(model) ???
# extract a solution a return it in the matrix form
solution = value.(x)
return Int64[ findfirst(isone, solution[i, j, :]) for i=1:9, j=1:9 ]
end
Objektivní funkce v prokládání exponenciálou.
@NLobjective(model, Min, sum((a[1] * exp(a[2] * ts[j]) + a[3] - ys[j])^2 for j = axes(ts, 1)))
Reference
Hlavní stránku balíčku JuMP naleznete zde. K prozkoumání doporučuji podívat se na obsáhlou dokumentaci.
Pro nadšené studenty FIT zajímající se o optimalizační problémy by mohla být relativně přístupná přehledná monografie Boyd, Vandenberghe: Convex Optimization, Cambridge University Press, 2004.