FIT ČVUT Course Pages
BI-JUL — Programování v jazyku Julia
Jdi na navigaci předmětu

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 2024/2025 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.

versioninfo()

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 f:M→Rf: \mathcal{M} \to \mathbb{R}, kde množina M\mathcal{M} je popsána pomocí funkcí gj:D→Rg_j: \mathcal{D} \to \mathbb{R}, j∈m^j\in \hat m, a hk:D→Rh_k: \mathcal{D} \to \mathbb{R}, k∈n^k\in\hat n, následovně

M={x∈D∣gj(x)=0, j∈m^, hk(x)≤0, k∈n^}.\mathcal{M} = \{x \in \mathcal{D} \mid g_j(x) = 0, \, j\in\hat m,\, h_k(x) \leq 0, \, k\in\hat n \}.

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 D\mathcal{D} může být otevřená podmnožina Rd\mathbb{R}^d, nebo třeba podmnožina Zd\mathbb{Z}^d, či dokonce podmnožina {0,1}d\{0,1\}^d.

Množina M\mathcal{M} 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í x∗∈Mx_* \in \mathcal{M} splňující f(x∗)=max⁡Mff(x_*) = \max_\mathcal{M} f (případně min⁡Mf\min_\mathcal{M} f) se nazývá "optimální řešení" a f(x∗)f(x_*) "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

2. Základní formulace problému ("modelu") v JuMP.jl

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 f(x,y)=x2−xy−y2f(x, y) = x^2 -xy - y^2 na obdélníku ⟨−1,1⟩×⟨−1,1⟩\langle -1, 1 \rangle \times \langle -1, 1 \rangle. 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 0x73170bbcc1a0>

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é xx a yy. 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 (alias for GenericVariableRef{Float64})
2 * x
$ 2 x $
x ^ 2
$ x^2 $

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)
$$ \begin{aligned} \min\quad & x^2 - x\times y - y^2\\ \text{Subject to} \quad & y \geq -1\\ & y \leq 1\\ & x \geq -1\\ & x \leq 1\\ \end{aligned} $$
latex_formulation(model)
$$ \begin{aligned} \min\quad & x^2 - x\times y - y^2\\ \text{Subject to} \quad & y \geq -1\\ & y \leq 1\\ & x \geq -1\\ & x \leq 1\\ \end{aligned} $$

Případně:

model
A JuMP Model
├ solver: Ipopt
├ objective_sense: MIN_SENSE
│ └ objective_function_type: QuadExpr
├ num_variables: 2
├ num_constraints: 4
│ ├ AffExpr in MOI.GreaterThan{Float64}: 1
│ ├ AffExpr in MOI.LessThan{Float64}: 1
│ ├ VariableRef in MOI.GreaterThan{Float64}: 1
│ └ VariableRef in MOI.LessThan{Float64}: 1
└ Names registered in the model
  └ :x, :y

Řešení spustíme metodou optimize!:

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.14.16, running with linear solver MUMPS 5.7.3.

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  0.0000000e+00 0.00e+00 2.00e-07  -1.7 0.00e+00   0.0 1.00e+00 1.00e+00   0
   2  0.0000000e+00 0.00e+00 1.50e-09  -3.8 0.00e+00   0.4 1.00e+00 1.00e+00   0
   3  0.0000000e+00 0.00e+00 1.84e-11  -5.7 0.00e+00  -0.1 1.00e+00 1.00e+00   0
   4  0.0000000e+00 0.00e+00 2.51e-14  -8.6 0.00e+00   0.4 1.00e+00 1.00e+00T  0

Number of Iterations....: 4

                                   (scaled)                 (unscaled)
Objective...............:   0.0000000000000000e+00    0.0000000000000000e+00
Dual infeasibility......:   2.5059035640133008e-14    2.5059035640133008e-14
Constraint violation....:   0.0000000000000000e+00    0.0000000000000000e+00
Variable bound violation:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   2.5059035596800808e-09    2.5059035596800808e-09
Overall NLP error.......:   2.5059035596800808e-09    2.5059035596800808e-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 seconds in IPOPT                               = 0.006

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
  Result count       : 1
  Termination status : LOCALLY_SOLVED
  Message from the solver:
  "Solve_Succeeded"

* Candidate solution (result #1)
  Primal status      : FEASIBLE_POINT
  Dual status        : FEASIBLE_POINT
  Objective value    : 0.00000e+00
  Dual objective value : -5.01176e-09

* Work counters
  Solve time (sec)   : 7.41196e-03
  Barrier iterations : 4

K optimálnímu řešení se dostaneme pomocí metody value.

value(x)
0.0
value(y)
0.0
value.([x, y])
2-element Vector{Float64}:
 0.0
 0.0

A k optimální hodnotě pomocí objective_value:

objective_value(model)
0.0

Vidíme, že se Ipopt řešič nedokázal dostat ze stacionárního bodu (0,0)(0,0), 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.14.16, running with linear solver MUMPS 5.7.3.

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
Variable bound 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             = 16
Number of objective gradient evaluations             = 10
Number of equality constraint evaluations            = 0
Number of inequality constraint evaluations          = 16
Number of equality constraint Jacobian evaluations   = 0
Number of inequality constraint Jacobian evaluations = 1
Number of Lagrangian Hessian evaluations             = 1
Total seconds in IPOPT                               = 0.009

EXIT: Optimal Solution Found.

Optimální řešení:

value.([x, y])
2-element Vector{Float64}:
 0.500000002821662
 1.0000000089982806
objective_value(model)
-1.2500000224957013

Kdybychom se vychýlili na druhou stranu, tak dokonvergujeme do jiného lokálního minima:

set_start_value(y, -0.1)
optimize!(model)
This is Ipopt version 3.14.16, running with linear solver MUMPS 5.7.3.

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  6.5900000e-02 0.00e+00 5.60e-01  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  9.0101477e-03 0.00e+00 1.32e-01  -1.0 1.24e-01   0.0 9.60e-01 1.00e+00f  1
   2 -1.9639654e-02 0.00e+00 2.05e-01  -1.7 7.70e-02   0.4 1.00e+00 1.00e+00f  1
   3 -1.1977756e+00 0.00e+00 7.52e-01  -1.7 9.49e-01  -0.1 1.00e+00 8.78e-01f  1
   4 -1.2176372e+00 0.00e+00 1.42e+00  -1.7 1.76e-01    -  1.00e+00 1.78e-01f  2
   5 -1.2299598e+00 0.00e+00 2.00e-07  -1.7 1.33e-01    -  1.00e+00 1.00e+00f  1
   6 -1.2499964e+00 0.00e+00 4.51e-05  -3.8 1.58e-02    -  9.97e-01 9.99e-01f  1
   7 -1.2499982e+00 0.00e+00 1.84e-11  -5.7 7.67e-04    -  1.00e+00 1.00e+00f  1
   8 -1.2500000e+00 0.00e+00 2.51e-14  -8.6 2.81e-06    -  1.00e+00 1.00e+00f  1

Number of Iterations....: 8

                                   (scaled)                 (unscaled)
Objective...............:  -1.2500000224965735e+00   -1.2500000224965735e+00
Dual infeasibility......:   2.5059035640133008e-14    2.5059035640133008e-14
Constraint violation....:   0.0000000000000000e+00    0.0000000000000000e+00
Variable bound violation:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   2.5230731469757146e-09    2.5230731469757146e-09
Overall NLP error.......:   2.5230731469757146e-09    2.5230731469757146e-09


Number of objective function evaluations             = 14
Number of objective gradient evaluations             = 9
Number of equality constraint evaluations            = 0
Number of inequality constraint evaluations          = 14
Number of equality constraint Jacobian evaluations   = 0
Number of inequality constraint Jacobian evaluations = 1
Number of Lagrangian Hessian evaluations             = 1
Total seconds in IPOPT                               = 0.005

EXIT: Optimal Solution Found.

Optimální řešení:

value.([x,y])
2-element Vector{Float64}:
 -0.5000000028103928
 -1.0000000089986294
objective_value(model)
-1.2500000224965735

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

Gurobi Error 10009: License expired 2023-12-11

Stacktrace:
  [1] _check_ret
    @ ~/.julia/packages/Gurobi/fnju1/src/MOI_wrapper/MOI_wrapper.jl:442 [inlined]
  [2] Gurobi.Env(params::Nothing; started::Bool, output_flag::Int64, memory_limit::Nothing)
    @ Gurobi ~/.julia/packages/Gurobi/fnju1/src/MOI_wrapper/MOI_wrapper.jl:176
  [3] Env (repeats 2 times)
    @ ~/.julia/packages/Gurobi/fnju1/src/MOI_wrapper/MOI_wrapper.jl:143 [inlined]
  [4] Gurobi.Optimizer(env::Nothing; enable_interrupts::Bool)
    @ Gurobi ~/.julia/packages/Gurobi/fnju1/src/MOI_wrapper/MOI_wrapper.jl:372
  [5] Optimizer
    @ ~/.julia/packages/Gurobi/fnju1/src/MOI_wrapper/MOI_wrapper.jl:366 [inlined]
  [6] Gurobi.Optimizer()
    @ Gurobi ~/.julia/packages/Gurobi/fnju1/src/MOI_wrapper/MOI_wrapper.jl:366
  [7] _instantiate_and_check(optimizer_constructor::Any)
    @ MathOptInterface ~/.julia/packages/MathOptInterface/gLl4d/src/instantiate.jl:94
  [8] instantiate(optimizer_constructor::Any; with_bridge_type::Type{Float64}, with_cache_type::Nothing)
    @ MathOptInterface ~/.julia/packages/MathOptInterface/gLl4d/src/instantiate.jl:175
  [9] set_optimizer(model::Model, optimizer_constructor::Any; add_bridges::Bool)
    @ JuMP ~/.julia/packages/JuMP/i68GU/src/optimizer_interface.jl:493
 [10] set_optimizer
    @ ~/.julia/packages/JuMP/i68GU/src/optimizer_interface.jl:486 [inlined]
 [11] Model(optimizer_factory::Any; add_bridges::Bool)
    @ JuMP ~/.julia/packages/JuMP/i68GU/src/JuMP.jl:195
 [12] Model(optimizer_factory::Any)
    @ JuMP ~/.julia/packages/JuMP/i68GU/src/JuMP.jl:187
 [13] top-level scope
    @ In[38]:3

Opět získáváme jedno z možných optimálních řešení:

value.([x,y])
objective_value(model)

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 SS. 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ě o1o_1, o2o_2 a o3o_3.
  • Péče o jednotku plochy oseté danou plodinou ho bude stát p1p_1, p2p_2 a p3p_3.
  • Výnos z jednotky plochy oseté danou plodinou bude v1v_1, v2v_2 a v3v_3.

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 x1x_1, x2x_2 a x3x_3. Jistě musí platit následující čtyři nerovnosti:

x1≥0, x2≥0, x3≥0, x1+x2+x3≤S.x_1 \geq 0, \ x_2 \geq 0, \ x_3 \geq 0, \ x_1 + x_2 + x_3 \leq S.

Snažíme se maximalizovat objektivní funkci (profit; k zjednodušení zápisu jsme přešli k vektorové notaci)

F(x1,x2,x3)=v⋅x−p⋅x−o⋅x=(v−p−o)⋅x.F(x_1, x_2, x_3) = v \cdot x - p \cdot x - o \cdot x = (v - p - o) \cdot x.

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,   0.9, 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)
$$ \begin{aligned} \max\quad & 3.5 x_{1} + 3.6 x_{2} + 3 x_{3}\\ \text{Subject to} \quad & x_{1} + x_{2} + x_{3} \leq 10\\ & x_{1} \geq 0\\ & x_{2} \geq 0\\ & x_{3} \geq 0\\ \end{aligned} $$

Spusťme řešení!

optimize!(model)
This is Ipopt version 3.14.16, running with linear solver MUMPS 5.7.3.

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  1.0099990e-01 0.00e+00 2.52e+00  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  1.4480398e+00 0.00e+00 3.56e+00  -1.0 3.99e-01    -  2.18e-01 1.00e+00f  1
   2  3.5160296e+01 0.00e+00 2.78e+00  -1.0 3.19e+01    -  1.03e-02 2.97e-01f  1
   3  3.5560614e+01 0.00e+00 2.35e+00  -1.0 8.39e+00    -  1.00e+00 1.57e-01f  2
   4  3.5711258e+01 0.00e+00 1.00e-06  -1.0 3.67e-01    -  1.00e+00 1.00e+00f  1
   5  3.5979796e+01 0.00e+00 9.21e-04  -2.5 9.66e-01    -  9.23e-01 1.00e+00f  1
   6  3.5999547e+01 0.00e+00 1.50e-09  -3.8 1.30e-01    -  1.00e+00 1.00e+00f  1
   7  3.5999995e+01 0.00e+00 1.84e-11  -5.7 1.79e-03    -  1.00e+00 1.00e+00f  1
   8  3.6000000e+01 0.00e+00 2.49e-14  -8.6 2.22e-05    -  1.00e+00 1.00e+00f  1

Number of Iterations....: 8

                                   (scaled)                 (unscaled)
Objective...............:  -3.6000000359478250e+01    3.6000000359478250e+01
Dual infeasibility......:   2.4900169555147560e-14    2.4900169555147560e-14
Constraint violation....:   0.0000000000000000e+00    0.0000000000000000e+00
Variable bound violation:   5.8225627789055691e-09    5.8225627789055691e-09
Complementarity.........:   2.5092896506099746e-09    2.5092896506099746e-09
Overall NLP error.......:   2.5092896506099746e-09    2.5092896506099746e-09


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 seconds in IPOPT                               = 0.004

EXIT: Optimal Solution Found.

Optimální řešení a hodnota:

value.(x)
3-element Vector{Float64}:
  1.5092896443320668e-8
 10.000000090033556
 -5.822562778905569e-9
objective_value(model)
36.00000035947825

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 "3535".

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 L>0L > 0, který je zavěšen v bodech (0,a)(0, a) a (D,a)(D, a), D,a>0D, a > 0, v homogenním gravitačním poli (v záporném směru osy yy). 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 nn hmotných bodů spojených nehmotným segmentem. Zajímá nás samozřejmě limita n→∞n \to \infty, 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)
$$ \begin{aligned} \min\quad & y_{1} + y_{2} + y_{3} + y_{4} + y_{5} + y_{6} + y_{7} + y_{8} + y_{9} + y_{10}\\ \text{Subject to} \quad & x_{1} = 0\\ & x_{10} = 5\\ & y_{1} = 2\\ & y_{10} = 2\\ & x_{1}^2 - 2 x_{2}\times x_{1} + x_{2}^2 + y_{1}^2 - 2 y_{2}\times y_{1} + y_{2}^2 = 0.6400000000000001\\ & x_{2}^2 - 2 x_{3}\times x_{2} + x_{3}^2 + y_{2}^2 - 2 y_{3}\times y_{2} + y_{3}^2 = 0.6400000000000001\\ & x_{3}^2 - 2 x_{4}\times x_{3} + x_{4}^2 + y_{3}^2 - 2 y_{4}\times y_{3} + y_{4}^2 = 0.6400000000000001\\ & x_{4}^2 - 2 x_{5}\times x_{4} + x_{5}^2 + y_{4}^2 - 2 y_{5}\times y_{4} + y_{5}^2 = 0.6400000000000001\\ & x_{5}^2 - 2 x_{6}\times x_{5} + x_{6}^2 + y_{5}^2 - 2 y_{6}\times y_{5} + y_{6}^2 = 0.6400000000000001\\ & x_{6}^2 - 2 x_{7}\times x_{6} + x_{7}^2 + y_{6}^2 - 2 y_{7}\times y_{6} + y_{7}^2 = 0.6400000000000001\\ & x_{7}^2 - 2 x_{8}\times x_{7} + x_{8}^2 + y_{7}^2 - 2 y_{8}\times y_{7} + y_{8}^2 = 0.6400000000000001\\ & x_{8}^2 - 2 x_{9}\times x_{8} + x_{9}^2 + y_{8}^2 - 2 y_{9}\times y_{8} + y_{9}^2 = 0.6400000000000001\\ & x_{9}^2 - 2 x_{10}\times x_{9} + x_{10}^2 + y_{9}^2 - 2 y_{10}\times y_{9} + y_{10}^2 = 0.6400000000000001\\ \end{aligned} $$

Pokusme se nalézt řešení.

optimize!(model)
This is Ipopt version 3.14.16, running with linear solver MUMPS 5.7.3.

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.8685030e+01 1.54e+00 3.29e+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.9976826e+01 1.45e+00 5.88e+07  -1.0 1.97e-01    -  1.00e+00 1.00e+00H  1
  11  1.9999793e+01 1.43e+00 1.11e+07  -1.0 2.61e-02    -  1.00e+00 1.00e+00h  1
  12  1.9999999e+01 1.39e+00 2.58e+06  -1.0 3.29e-02    -  1.00e+00 1.00e+00h  1
  13  2.0000001e+01 1.36e+00 2.51e+06  -1.0 3.53e-02    -  1.00e+00 1.00e+00h  1
  14  2.0000001e+01 1.32e+00 2.48e+06  -1.0 3.73e-02    -  1.00e+00 1.00e+00h  1
  15  2.0000001e+01 1.28e+00 2.44e+06  -1.0 3.82e-02    -  1.00e+00 1.00e+00h  1
  16  2.0000001e+01 1.25e+00 2.41e+06  -1.0 3.87e-02    -  1.00e+00 1.00e+00h  1
  17  2.0000001e+01 1.21e+00 2.37e+06  -1.0 3.88e-02    -  1.00e+00 1.00e+00h  1
  18  2.0000001e+01 1.17e+00 2.33e+06  -1.0 3.88e-02    -  1.00e+00 1.00e+00h  1
  19  2.0000001e+01 1.13e+00 2.29e+06  -1.0 3.86e-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.0000001e+01 1.09e+00 2.25e+06  -1.0 3.83e-02    -  1.00e+00 1.00e+00h  1
  21  2.0000001e+01 1.05e+00 2.20e+06  -1.0 3.79e-02    -  1.00e+00 1.00e+00h  1
  22  2.0000001e+01 1.02e+00 2.16e+06  -1.0 3.75e-02    -  1.00e+00 1.00e+00h  1
  23  2.0000001e+01 9.78e-01 2.11e+06  -1.0 3.71e-02    -  1.00e+00 1.00e+00h  1
  24  2.0000001e+01 9.25e-01 4.52e+06  -1.7 5.36e-02    -  1.00e+00 1.00e+00h  1
  25  2.0000000e+01 8.72e-01 4.40e+06  -1.7 5.30e-02    -  1.00e+00 1.00e+00h  1
  26  2.0000000e+01 8.20e-01 4.25e+06  -1.7 5.18e-02    -  1.00e+00 1.00e+00h  1
  27  2.0000000e+01 7.69e-01 4.09e+06  -1.7 5.05e-02    -  1.00e+00 1.00e+00h  1
  28  2.0000000e+01 7.20e-01 3.94e+06  -1.7 4.91e-02    -  1.00e+00 1.00e+00h  1
  29  2.0000000e+01 6.73e-01 3.78e+06  -1.7 4.77e-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.0000000e+01 6.26e-01 3.62e+06  -1.7 4.62e-02    -  1.00e+00 1.00e+00h  1
  31  2.0000000e+01 5.98e-01 3.47e+06  -1.7 4.47e-02    -  1.00e+00 1.00e+00h  1
  32  2.0000000e+01 5.96e-01 3.32e+06  -1.7 4.31e-02    -  1.00e+00 1.00e+00h  1
  33  2.0000000e+01 5.93e-01 3.17e+06  -1.7 4.16e-02    -  1.00e+00 1.00e+00h  1
  34  2.0000000e+01 5.91e-01 3.02e+06  -1.7 4.00e-02    -  1.00e+00 1.00e+00h  1
  35  2.0000000e+01 5.88e-01 2.87e+06  -1.7 3.84e-02    -  1.00e+00 1.00e+00h  1
  36  2.0000000e+01 5.86e-01 2.73e+06  -1.7 3.67e-02    -  1.00e+00 1.00e+00h  1
  37  2.0000000e+01 5.83e-01 2.60e+06  -1.7 3.51e-02    -  1.00e+00 1.00e+00h  1
  38  2.0000000e+01 5.81e-01 2.47e+06  -1.7 3.35e-02    -  1.00e+00 1.00e+00h  1
  39  2.0000000e+01 5.79e-01 2.34e+06  -1.7 3.19e-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.0000000e+01 5.77e-01 2.21e+06  -1.7 3.04e-02    -  1.00e+00 1.00e+00h  1
  41  2.0000000e+01 5.74e-01 2.10e+06  -1.7 2.89e-02    -  1.00e+00 1.00e+00h  1
  42  2.0000000e+01 5.72e-01 1.98e+06  -1.7 2.74e-02    -  1.00e+00 1.00e+00h  1
  43  2.0000000e+01 5.70e-01 1.87e+06  -1.7 2.60e-02    -  1.00e+00 1.00e+00h  1
  44  2.0000000e+01 5.68e-01 1.77e+06  -1.7 2.46e-02    -  1.00e+00 1.00e+00h  1
  45  2.0000000e+01 5.66e-01 1.67e+06  -1.7 2.32e-02    -  1.00e+00 1.00e+00h  1
  46  2.0000000e+01 5.64e-01 1.57e+06  -1.7 2.20e-02    -  1.00e+00 1.00e+00h  1
  47  2.0000000e+01 5.62e-01 1.48e+06  -1.7 2.07e-02    -  1.00e+00 1.00e+00h  1
  48  2.0000000e+01 5.60e-01 1.40e+06  -1.7 1.96e-02    -  1.00e+00 1.00e+00h  1
  49  2.0000000e+01 5.58e-01 1.32e+06  -1.7 1.85e-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.0000000e+01 5.57e-01 1.24e+06  -1.7 1.74e-02    -  1.00e+00 1.00e+00h  1
  51  2.0000000e+01 5.55e-01 1.16e+06  -1.7 1.64e-02    -  1.00e+00 1.00e+00h  1
  52  2.0000000e+01 5.53e-01 1.10e+06  -1.7 1.55e-02    -  1.00e+00 1.00e+00h  1
  53  2.0000000e+01 5.51e-01 1.08e+06  -1.7 1.46e-02    -  1.00e+00 1.00e+00h  1
  54  2.0000000e+01 5.50e-01 1.05e+06  -1.7 1.37e-02    -  1.00e+00 1.00e+00h  1
  55  2.0000000e+01 5.48e-01 1.03e+06  -1.7 1.29e-02    -  1.00e+00 1.00e+00h  1
  56  2.0000000e+01 5.47e-01 9.99e+05  -1.7 1.22e-02    -  1.00e+00 1.00e+00h  1
  57  2.0000000e+01 5.45e-01 9.72e+05  -1.7 1.15e-02    -  1.00e+00 1.00e+00h  1
  58  2.0000000e+01 5.43e-01 9.44e+05  -1.7 1.08e-02    -  1.00e+00 1.00e+00h  1
  59  2.0000000e+01 5.42e-01 9.17e+05  -1.7 1.02e-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.0000000e+01 5.40e-01 8.90e+05  -1.7 9.59e-03    -  1.00e+00 1.00e+00h  1
  61  2.0000000e+01 5.39e-01 8.63e+05  -1.7 9.04e-03    -  1.00e+00 1.00e+00h  1
  62  2.0000000e+01 5.37e-01 8.37e+05  -1.7 8.52e-03    -  1.00e+00 1.00e+00h  1
  63  2.0000000e+01 5.36e-01 8.12e+05  -1.7 8.03e-03    -  1.00e+00 1.00e+00h  1
  64  2.0000000e+01 5.34e-01 7.87e+05  -1.7 7.58e-03    -  1.00e+00 1.00e+00h  1
  65  2.0000000e+01 5.33e-01 7.64e+05  -1.7 7.15e-03    -  1.00e+00 1.00e+00h  1
  66  2.0000000e+01 5.32e-01 7.41e+05  -1.7 6.75e-03    -  1.00e+00 1.00e+00h  1
  67  2.0000000e+01 5.31e-01 5.50e+05  -1.7 6.37e-03    -  1.00e+00 5.00e-01h  2
  68r 2.0000000e+01 5.31e-01 9.99e+02  -0.3 0.00e+00    -  0.00e+00 4.77e-07R 22
  69r 2.0000001e+01 6.19e-01 9.94e+02  -0.3 8.13e+01    -  6.47e-03 5.11e-03f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  70r 2.0000001e+01 6.08e-01 9.91e+02  -0.3 5.06e+00   2.0 1.00e+00 4.27e-03f  1
  71r 2.0000001e+01 5.81e-01 8.68e+02  -0.3 8.22e-01   2.4 1.00e+00 1.60e-01f  1
  72r 1.9999982e+01 6.06e-01 4.46e+02  -0.3 3.59e-01   2.9 1.00e+00 6.88e-01f  1
  73r 1.9999981e+01 6.10e-01 1.15e+03  -0.3 8.67e-02   4.2 1.00e+00 8.09e-01f  1
  74r 1.9999961e+01 6.13e-01 4.07e+02  -0.3 3.27e-02   3.7 1.00e+00 7.54e-01f  1
  75r 1.9999946e+01 6.15e-01 1.42e+02  -0.3 1.05e-02   4.1 1.00e+00 1.00e+00f  1
  76r 1.9999748e+01 6.21e-01 9.71e+01  -0.3 2.11e-02   3.7 1.00e+00 1.00e+00f  1
  77r 1.9999653e+01 6.23e-01 1.04e+02  -0.3 8.71e-03   4.1 1.00e+00 1.00e+00f  1
  78r 1.9995991e+01 6.22e-01 1.69e+02  -0.3 4.14e-02   3.6 1.00e+00 1.00e+00f  1
  79r 1.9994578e+01 6.16e-01 2.13e+02  -0.3 2.10e-02   4.0 1.00e+00 1.00e+00f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  80r 1.9974084e+01 3.21e-01 6.87e+02  -0.3 3.39e-01   3.6 1.00e+00 1.00e+00f  1
  81  1.9962279e+01 1.09e-01 2.77e+07  -1.7 8.25e-02   8.3 1.00e+00 1.00e+00h  1
  82  1.9953554e+01 8.16e-02 1.65e+07  -1.7 2.88e-02   8.7 1.00e+00 1.00e+00h  1
  83  1.9927126e+01 9.51e-02 1.32e+07  -1.7 7.26e-02   8.3 1.00e+00 1.00e+00h  1
  84  1.9919068e+01 7.82e-02 1.60e+07  -1.7 3.31e-02   8.7 1.00e+00 1.00e+00h  1
  85  1.9902171e+01 2.33e-02 1.45e+07  -1.7 8.97e-02   8.2 1.00e+00 1.00e+00h  1
  86  1.9897614e+01 1.28e-02 1.30e+07  -1.7 3.01e-02   8.6 1.00e+00 1.00e+00h  1
  87  1.9836231e+01 1.05e-01 2.01e+07  -1.7 1.58e-01   8.2 1.00e+00 1.00e+00h  1
  88  1.9799203e+01 6.53e-02 1.78e+07  -1.7 5.00e-02   8.6 1.00e+00 1.00e+00h  1
  89  1.9690072e+01 1.14e-01 1.03e+07  -1.7 9.44e-02   8.1 1.00e+00 1.00e+00h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  90  1.9624383e+01 1.12e-01 2.30e+06  -1.7 4.11e-02   7.6 1.00e+00 1.00e+00h  1
  91  1.9803196e+01 1.03e-01 1.61e+06  -1.7 1.59e-01   7.2 1.00e+00 1.00e+00h  1
  92  1.9866256e+01 2.15e-02 1.94e+06  -1.7 6.12e-02   7.6 1.00e+00 1.00e+00h  1
  93  1.9882079e+01 6.97e-03 2.02e+06  -1.7 2.03e-02   8.0 1.00e+00 1.00e+00h  1
  94  1.9888607e+01 4.59e-03 3.59e+05  -1.7 1.66e-01   7.5 1.00e+00 1.00e+00H  1
  95  1.9894974e+01 4.93e-03 1.26e+05  -1.7 1.13e-02   7.0 1.00e+00 1.00e+00h  1
  96  1.9893199e+01 2.47e-03 1.13e+04  -1.7 3.06e-03   6.6 1.00e+00 1.00e+00h  1
  97  1.9889481e+01 7.21e-05 6.75e+03  -1.7 5.46e-03   6.1 1.00e+00 1.00e+00h  1
  98  1.9889320e+01 1.37e-04 8.79e+01  -1.7 2.17e-04   5.6 1.00e+00 1.00e+00h  1
  99  1.9889453e+01 2.44e-06 3.70e+01  -1.7 2.68e-04   5.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.9889291e+01 7.77e-07 1.17e+00  -1.7 2.54e-05   4.7 1.00e+00 1.00e+00h  1
 101  1.9888791e+01 8.43e-09 1.03e+00  -1.7 6.72e-05   4.2 1.00e+00 1.00e+00h  1
 102  1.9887294e+01 4.08e-08 1.04e+00  -1.7 2.03e-04   3.7 1.00e+00 1.00e+00f  1
 103  1.9882803e+01 3.70e-07 1.04e+00  -1.7 6.10e-04   3.2 1.00e+00 1.00e+00f  1
 104  1.9869322e+01 3.34e-06 1.04e+00  -1.7 1.83e-03   2.8 1.00e+00 1.00e+00f  1
 105  1.9828829e+01 3.02e-05 1.04e+00  -1.7 5.49e-03   2.3 1.00e+00 1.00e+00f  1
 106  1.9706903e+01 2.75e-04 1.04e+00  -1.7 1.65e-02   1.8 1.00e+00 1.00e+00f  1
 107  1.9337358e+01 2.57e-03 1.06e+00  -1.7 5.04e-02   1.3 1.00e+00 1.00e+00f  1
 108  1.8202727e+01 2.37e-02 1.03e+00  -1.7 1.54e-01   0.8 1.00e+00 1.00e+00f  1
 109  1.4758607e+01 2.13e-01 1.16e+00  -1.7 4.75e-01   0.4 1.00e+00 1.00e+00f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
 110  3.7168507e+00 3.35e+00 4.10e+00  -1.7 2.09e+00  -0.1 1.00e+00 1.00e+00f  1
 111  5.5254855e+00 7.47e-01 1.22e+01  -1.7 8.38e-01   1.2 1.00e+00 1.00e+00h  1
 112  6.6807345e+00 2.67e-01 1.31e+01  -1.7 4.69e-01   1.6 1.00e+00 1.00e+00h  1
 113  7.8533879e+00 2.14e-02 2.79e+00  -1.7 2.41e-01    -  1.00e+00 1.00e+00h  1
 114  7.1452945e+00 4.82e-02 5.16e-01  -1.7 3.33e-01    -  1.00e+00 1.00e+00f  1
 115  7.2255909e+00 3.28e-02 1.05e-01  -1.7 1.37e-01    -  1.00e+00 1.00e+00h  1
 116  7.3670454e+00 2.07e-03 3.04e-03  -2.5 5.64e-02    -  1.00e+00 1.00e+00h  1
 117  7.3730869e+00 1.90e-06 1.40e-05  -3.8 1.35e-03    -  1.00e+00 1.00e+00h  1
 118  7.3730955e+00 4.85e-12 5.15e-11  -8.6 2.86e-06    -  1.00e+00 1.00e+00h  1

Number of Iterations....: 118

                                   (scaled)                 (unscaled)
Objective...............:   7.3730955000157827e+00    7.3730955000157827e+00
Dual infeasibility......:   5.1451620741715942e-11    5.1451620741715942e-11
Constraint violation....:   4.8540060859636469e-12    4.8540060859636469e-12
Variable bound violation:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   0.0000000000000000e+00    0.0000000000000000e+00
Overall NLP error.......:   5.1451620741715942e-11    5.1451620741715942e-11


Number of objective function evaluations             = 148
Number of objective gradient evaluations             = 108
Number of equality constraint evaluations            = 148
Number of inequality constraint evaluations          = 0
Number of equality constraint Jacobian evaluations   = 120
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations             = 118
Total seconds in IPOPT                               = 0.052

EXIT: Optimal Solution Found.

solution_summary(model)
* Solver : Ipopt

* Status
  Result count       : 1
  Termination status : LOCALLY_SOLVED
  Message from the solver:
  "Solve_Succeeded"

* Candidate solution (result #1)
  Primal status      : FEASIBLE_POINT
  Dual status        : FEASIBLE_POINT
  Objective value    : 7.37310e+00
  Dual objective value : 1.87720e+01

* Work counters
  Solve time (sec)   : 5.31480e-02
  Barrier iterations : 118

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 0x731702b62f00>

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 11 až 99 v 9×99 \times 9 matici (Mi,j)i,j=19(M_{i,j})_{i,j=1}^9 si představme jedničky a nuly v 9×9×99\times9\times9 krychli (Ki,j,k)i,j,k=19(K_{i,j,k})_{i,j,k=1}^9. Krychle je vyplněna nulami, jen a pouze je-li v původní matici hodnota Mi,jM_{i,j}, pak tomu odpovídá hodnota Ki,j,Mi,j=1K_{i,j,M_{i,j}} = 1.

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 3×33\times3 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
│   Result count       : 1
│   Termination status : OPTIMAL
│   Message from the solver:
│   "Solution is optimal"
│ 
│ * Candidate solution (result #1)
│   Primal status      : FEASIBLE_POINT
│   Dual status        : NO_SOLUTION
│   Objective value    : 0.00000e+00
│   Objective bound    : 0.00000e+00
│   Relative gap       : 0.00000e+00
│ 
│ * Work counters
└   Solve time (sec)   : 4.49839e-02

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
multi = zeros(Int64, 9, 9)
9×9 Matrix{Int64}:
 0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0
sudoku_solver(multi)
┌ Info: * Solver : GLPK
│ 
│ * Status
│   Result count       : 1
│   Termination status : OPTIMAL
│   Message from the solver:
│   "Solution is optimal"
│ 
│ * Candidate solution (result #1)
│   Primal status      : FEASIBLE_POINT
│   Dual status        : NO_SOLUTION
│   Objective value    : 0.00000e+00
│   Objective bound    : 0.00000e+00
│   Relative gap       : 0.00000e+00
│ 
│ * Work counters
└   Solve time (sec)   : 4.72606e-01

9×9 Matrix{Int64}:
 5  2  7  6  4  1  8  3  9
 1  9  3  5  7  8  4  2  6
 6  8  4  2  9  3  7  1  5
 4  5  8  3  1  7  9  6  2
 7  6  1  8  2  9  3  5  4
 9  3  2  4  6  5  1  8  7
 2  7  6  1  3  4  5  9  8
 8  1  9  7  5  6  2  4  3
 3  4  5  9  8  2  6  7  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 (yi)i=1n(y_i)_{i=1}^n nějaké veličiny yy v časových okamžicích (ti)i=1n(t_i)_{i=1}^n (nezávisle proměnná nutně nemusí být čas tt). Dále máme podezření, že hypotetická závislost y(t)y(t) se chová exponenciálně, obecně třeba ve tvaru αeβt+γ\alpha e^{\beta t} + \gamma, kde α,β,γ\alpha,\beta,\gamma 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 F(α,β,γ)=∑i=1n(αeβti+γ−yi)2,α,β,γ∈R,F(\alpha, \beta, \gamma) = \sum_{i=1}^n \Big( \alpha e^{\beta t_i} + \gamma - y_i \Big)^2, \quad \alpha,\beta,\gamma \in \mathbb{R}, 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)

    if optimizer == Ipopt.Optimizer
        set_optimizer_attribute(model, "max_iter", max_iter)
    end
    
    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 is Ipopt version 3.14.16, running with linear solver MUMPS 5.7.3.

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.5576370e+04 0.00e+00 1.00e+02  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  9.5536297e+04 0.00e+00 1.00e+02  -1.0 1.00e-02   4.0 1.00e+00 1.00e+00f  1
   2  9.5414395e+04 0.00e+00 1.01e+02  -1.0 3.04e-02   3.5 1.00e+00 1.00e+00f  1
   3  9.4959057e+04 0.00e+00 1.21e+02  -1.0 1.07e-01   3.0 1.00e+00 1.00e+00f  1
   4  9.4727161e+04 0.00e+00 1.35e+02  -1.0 4.54e-02   3.5 1.00e+00 1.00e+00f  1
   5  8.9075180e+04 0.00e+00 1.07e+03  -1.0 4.00e-01   3.0 1.00e+00 1.00e+00f  1
   6  8.7052117e+04 0.00e+00 1.42e+03  -1.0 6.52e-02   4.3 1.00e+00 1.00e+00f  1
   7  9.1146330e+03 0.00e+00 6.95e+03  -1.0 1.31e+00   3.8 1.00e+00 5.00e-01f  2
   8  2.0624975e+03 0.00e+00 5.09e+03  -1.0 1.43e+01    -  1.00e+00 5.00e-01f  2
   9  6.1547208e+02 0.00e+00 3.16e+02  -1.0 9.15e-02   3.4 1.00e+00 1.00e+00f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  10  5.8673084e+02 0.00e+00 5.73e+01  -1.0 3.64e-02   2.9 1.00e+00 1.00e+00f  1
  11  4.1450956e+02 0.00e+00 1.07e+03  -1.0 5.08e+00    -  1.00e+00 2.50e-01f  3
  12  1.7768208e+02 0.00e+00 1.65e+02  -1.0 1.91e+00    -  1.00e+00 1.00e+00f  1
  13  1.3409469e+02 0.00e+00 6.30e+02  -1.0 9.85e-01    -  1.00e+00 5.00e-01f  2
  14  8.5824125e+01 0.00e+00 9.06e+01  -1.0 2.35e-01    -  1.00e+00 1.00e+00f  1
  15  7.1830648e+01 0.00e+00 1.95e+02  -1.0 4.94e-01    -  1.00e+00 5.00e-01f  2
  16  6.2720633e+01 0.00e+00 9.45e+01  -1.0 2.02e-01    -  1.00e+00 1.00e+00f  1
  17  6.0328974e+01 0.00e+00 3.82e+01  -1.0 1.23e-01    -  1.00e+00 1.00e+00f  1
  18  6.0033084e+01 0.00e+00 5.58e+00  -1.0 4.68e-02    -  1.00e+00 1.00e+00f  1
  19  6.0025651e+01 0.00e+00 1.99e-01  -1.0 8.74e-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.0025644e+01 0.00e+00 1.55e-04  -2.5 2.46e-04    -  1.00e+00 1.00e+00f  1
  21  6.0025644e+01 0.00e+00 1.79e-10  -5.7 2.53e-07    -  1.00e+00 1.00e+00f  1

Number of Iterations....: 21

                                   (scaled)                 (unscaled)
Objective...............:   2.9953186524458988e+00    6.0025644095341846e+01
Dual infeasibility......:   1.7931934547028065e-10    3.5935272535425611e-09
Constraint violation....:   0.0000000000000000e+00    0.0000000000000000e+00
Variable bound violation:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   0.0000000000000000e+00    0.0000000000000000e+00
Overall NLP error.......:   1.7931934547028065e-10    3.5935272535425611e-09


Number of objective function evaluations             = 48
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 seconds in IPOPT                               = 0.010

EXIT: Optimal Solution Found.

┌ Info: * Solver : Ipopt
│ 
│ * Status
│   Result count       : 1
│   Termination status : LOCALLY_SOLVED
│   Message from the solver:
│   "Solve_Succeeded"
│ 
│ * Candidate solution (result #1)
│   Primal status      : FEASIBLE_POINT
│   Dual status        : FEASIBLE_POINT
│   Objective value    : 6.00256e+01
│   Dual objective value : 0.00000e+00
│ 
│ * Work counters
│   Solve time (sec)   : 1.06499e-02
└   Barrier iterations : 21

(60.025644095341846, [0.30267930857692293, 1.577153411335236, -2.321617556606973])

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.14.16, running with linear solver MUMPS 5.7.3.

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.4536867e+03 0.00e+00 1.00e+02  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  6.4216977e+03 0.00e+00 9.99e+01  -1.0 9.99e-03   4.0 1.00e+00 1.00e+00f  1
   2  6.3251513e+03 0.00e+00 1.01e+02  -1.0 3.02e-02   3.5 1.00e+00 1.00e+00f  1
   3  5.9738009e+03 0.00e+00 1.18e+02  -1.0 1.05e-01   3.0 1.00e+00 1.00e+00f  1
   4  5.7937940e+03 0.00e+00 1.32e+02  -1.0 4.43e-02   3.5 1.00e+00 1.00e+00f  1
   5  3.3614545e+03 0.00e+00 5.62e+02  -1.0 6.12e-01   3.0 1.00e+00 5.00e-01f  2
   6  2.6384275e+03 0.00e+00 2.70e+02  -1.0 8.58e-01   2.5 1.00e+00 1.25e-01f  4
   7  2.0651704e+03 0.00e+00 5.02e+02  -1.0 1.15e+00   2.0 1.00e+00 5.00e-01f  2
   8  1.7241816e+01 0.00e+00 1.74e+00  -1.0 1.45e+01    -  1.00e+00 5.00e-01f  2
   9  1.2253538e+01 0.00e+00 9.83e+00  -1.0 9.55e-01    -  1.00e+00 5.00e-01f  2
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  10  1.0507711e+01 0.00e+00 1.40e+00  -1.0 1.65e-02   1.6 1.00e+00 1.00e+00f  1
  11  1.0304187e+01 0.00e+00 2.76e-01  -1.0 2.27e-02   1.1 1.00e+00 1.00e+00f  1
  12  1.0035460e+01 0.00e+00 6.28e+00  -1.7 6.81e-01    -  1.00e+00 5.00e-01f  2
  13  9.2646387e+00 0.00e+00 2.36e+00  -1.7 2.39e-01    -  1.00e+00 1.00e+00f  1
  14  9.1257629e+00 0.00e+00 1.03e-01  -1.7 2.76e-02   0.6 1.00e+00 1.00e+00f  1
  15  9.0350066e+00 0.00e+00 7.33e+00  -2.5 8.27e-01    -  1.00e+00 5.00e-01f  2
  16  8.4976265e+00 0.00e+00 9.82e-01  -2.5 2.20e-01    -  1.00e+00 1.00e+00f  1
  17  8.3004208e+00 0.00e+00 1.73e+00  -2.5 4.75e+00    -  1.00e+00 6.25e-02f  5
  18  8.1370863e+00 0.00e+00 5.87e+00  -2.5 4.91e-01    -  1.00e+00 1.00e+00f  1
  19  7.9022141e+00 0.00e+00 3.87e-01  -2.5 3.08e-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  7.8178613e+00 0.00e+00 7.00e+00  -2.5 1.22e+00    -  1.00e+00 5.00e-01f  2
  21  7.6327141e+00 0.00e+00 2.92e-01  -2.5 3.53e-01    -  1.00e+00 1.00e+00f  1
  22  7.6059859e+00 0.00e+00 9.97e+00  -2.5 1.68e+00    -  1.00e+00 5.00e-01f  2
  23  7.4251453e+00 0.00e+00 2.61e-01  -2.5 3.78e-01    -  1.00e+00 1.00e+00f  1
  24  7.3726306e+00 0.00e+00 7.39e+00  -2.5 3.41e+00    -  1.00e+00 2.50e-01f  3
  25  7.2784075e+00 0.00e+00 8.69e-01  -2.5 5.10e-01    -  1.00e+00 1.00e+00f  1
  26  7.2409794e+00 0.00e+00 8.12e+00  -2.5 1.97e+00    -  1.00e+00 5.00e-01f  2
  27  7.1674969e+00 0.00e+00 9.00e-01  -2.5 5.88e-01    -  1.00e+00 1.00e+00f  1
  28  7.1420082e+00 0.00e+00 9.53e+00  -2.5 2.39e+00    -  1.00e+00 5.00e-01f  2
  29  7.0807127e+00 0.00e+00 7.26e-01  -2.5 6.38e-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  7.0748815e+00 0.00e+00 1.36e+01  -2.5 3.17e+00    -  1.00e+00 5.00e-01f  2
  31  7.0097032e+00 0.00e+00 4.39e-02  -2.5 6.29e-01    -  1.00e+00 1.00e+00f  1
  32  6.9911690e+00 0.00e+00 9.39e+00  -2.5 5.96e+00    -  1.00e+00 2.50e-01f  3
  33  6.9572991e+00 0.00e+00 1.49e+00  -2.5 8.83e-01    -  1.00e+00 1.00e+00f  1
  34  6.9429392e+00 0.00e+00 1.01e+01  -2.5 3.30e+00    -  1.00e+00 5.00e-01f  2
  35  6.9154262e+00 0.00e+00 1.56e+00  -2.5 9.83e-01    -  1.00e+00 1.00e+00f  1
  36  6.9043274e+00 0.00e+00 1.11e+01  -2.5 3.78e+00    -  1.00e+00 5.00e-01f  2
  37  6.8812536e+00 0.00e+00 1.50e+00  -2.5 1.06e+00    -  1.00e+00 1.00e+00f  1
  38  6.8739017e+00 0.00e+00 1.33e+01  -2.5 4.51e+00    -  1.00e+00 5.00e-01f  2
  39  6.8527895e+00 0.00e+00 1.10e+00  -2.5 1.08e+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  6.8440160e+00 0.00e+00 5.69e+00  -2.5 6.16e+00    -  1.00e+00 2.50e-01f  3
  41  6.8326552e+00 0.00e+00 8.79e+00  -2.5 2.24e+00    -  1.00e+00 1.00e+00f  1
  42  6.8206487e+00 0.00e+00 4.10e+00  -2.5 1.73e+00    -  1.00e+00 1.00e+00f  1
  43  6.8172678e+00 0.00e+00 1.81e+01  -2.5 3.32e+00    -  1.00e+00 1.00e+00f  1
  44  6.8009942e+00 0.00e+00 6.75e-01  -2.5 1.21e+00    -  1.00e+00 1.00e+00f  1
  45  6.7954442e+00 0.00e+00 1.10e+01  -2.5 1.07e+01    -  1.00e+00 2.50e-01f  3
  46  6.7866712e+00 0.00e+00 3.47e+00  -2.5 1.89e+00    -  1.00e+00 1.00e+00f  1
  47  6.7817715e+00 0.00e+00 9.12e+00  -2.5 4.91e+00    -  1.00e+00 5.00e-01f  2
  48  6.7748319e+00 0.00e+00 5.88e+00  -2.5 2.45e+00    -  1.00e+00 1.00e+00f  1
  49  6.7702140e+00 0.00e+00 1.44e+01  -2.5 3.69e+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  6.7631968e+00 0.00e+00 2.54e+00  -2.5 1.95e+00    -  1.00e+00 1.00e+00f  1
  51  6.7605028e+00 0.00e+00 1.63e+01  -2.5 7.80e+00    -  1.00e+00 5.00e-01f  2
  52  6.7540706e+00 0.00e+00 2.17e+00  -2.5 1.97e+00    -  1.00e+00 1.00e+00f  1
  53  6.7531334e+00 0.00e+00 2.21e+01  -2.5 9.69e+00    -  1.00e+00 5.00e-01f  2
  54  6.7458686e+00 0.00e+00 1.07e+00  -2.5 1.77e+00    -  1.00e+00 1.00e+00f  1
  55  6.7433562e+00 0.00e+00 1.36e+01  -2.5 1.61e+01    -  1.00e+00 2.50e-01f  3
  56  6.7392672e+00 0.00e+00 4.31e+00  -2.5 2.77e+00    -  1.00e+00 1.00e+00f  1
  57  6.7369677e+00 0.00e+00 1.14e+01  -2.5 7.38e+00    -  1.00e+00 5.00e-01f  2
  58  6.7336330e+00 0.00e+00 6.87e+00  -2.5 3.49e+00    -  1.00e+00 1.00e+00f  1
  59  6.7316684e+00 0.00e+00 1.92e+01  -2.5 5.61e+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  6.7279581e+00 0.00e+00 2.52e+00  -2.5 2.57e+00    -  1.00e+00 1.00e+00f  1
  61  6.7277580e+00 0.00e+00 2.72e+01  -2.5 1.34e+01    -  1.00e+00 5.00e-01f  2
  62  6.7231593e+00 0.00e+00 1.08e+00  -2.5 2.19e+00    -  1.00e+00 1.00e+00f  1
  63  6.7219668e+00 0.00e+00 1.98e+01  -2.5 2.41e+01    -  1.00e+00 2.50e-01f  3
  64  6.7191868e+00 0.00e+00 2.98e+00  -2.5 2.99e+00    -  1.00e+00 1.00e+00f  1
  65  6.7186742e+00 0.00e+00 2.61e+01  -2.5 1.44e+01    -  1.00e+00 5.00e-01f  2
  66  6.7155359e+00 0.00e+00 1.66e+00  -2.5 2.65e+00    -  1.00e+00 1.00e+00f  1
  67  6.7143647e+00 0.00e+00 1.55e+01  -2.5 2.31e+01    -  1.00e+00 2.50e-01f  3
  68  6.7125243e+00 0.00e+00 6.25e+00  -2.5 4.32e+00    -  1.00e+00 1.00e+00f  1
  69  6.7114552e+00 0.00e+00 1.20e+01  -2.5 9.71e+00    -  1.00e+00 5.00e-01f  2
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  70  6.7099499e+00 0.00e+00 1.14e+01  -2.5 5.80e+00    -  1.00e+00 1.00e+00f  1
  71  6.7085744e+00 0.00e+00 1.32e+01  -2.5 6.33e+00    -  1.00e+00 1.00e+00f  1
  72  6.7072085e+00 0.00e+00 1.04e+01  -2.5 5.84e+00    -  1.00e+00 1.00e+00f  1
  73  6.7060884e+00 0.00e+00 1.72e+01  -2.5 7.47e+00    -  1.00e+00 1.00e+00f  1
  74  6.7047729e+00 0.00e+00 6.67e+00  -2.5 5.06e+00    -  1.00e+00 1.00e+00f  1
  75  6.7040072e+00 0.00e+00 1.36e+01  -2.5 1.19e+01    -  1.00e+00 5.00e-01f  2
  76  6.7029150e+00 0.00e+00 1.16e+01  -2.5 6.66e+00    -  1.00e+00 1.00e+00f  1
  77  6.7019653e+00 0.00e+00 1.64e+01  -2.5 7.92e+00    -  1.00e+00 1.00e+00f  1
  78  6.7009361e+00 0.00e+00 8.76e+00  -2.5 6.15e+00    -  1.00e+00 1.00e+00f  1
  79  6.7005168e+00 0.00e+00 3.04e+01  -2.5 1.10e+01    -  1.00e+00 1.00e+00f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  80  6.6991430e+00 0.00e+00 2.46e+00  -2.5 3.98e+00    -  1.00e+00 1.00e+00f  1
  81  6.6985884e+00 0.00e+00 1.71e+01  -2.5 3.21e+01    -  1.00e+00 2.50e-01f  3
  82  6.6977494e+00 0.00e+00 9.42e+00  -2.5 6.85e+00    -  1.00e+00 1.00e+00f  1
  83  6.6973717e+00 0.00e+00 3.09e+01  -2.5 1.20e+01    -  1.00e+00 1.00e+00f  1
  84  6.6962819e+00 0.00e+00 2.84e+00  -2.5 4.49e+00    -  1.00e+00 1.00e+00f  1
  85  6.6958237e+00 0.00e+00 1.66e+01  -2.5 3.37e+01    -  1.00e+00 2.50e-01f  3
  86  6.6951479e+00 0.00e+00 1.17e+01  -2.5 8.10e+00    -  1.00e+00 1.00e+00f  1
  87  6.6946477e+00 0.00e+00 2.40e+01  -2.5 1.14e+01    -  1.00e+00 1.00e+00f  1
  88  6.6939440e+00 0.00e+00 5.94e+00  -2.5 6.32e+00    -  1.00e+00 1.00e+00f  1
  89  6.6935966e+00 0.00e+00 2.42e+01  -2.5 2.19e+01    -  1.00e+00 5.00e-01f  2
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  90  6.6929565e+00 0.00e+00 6.27e+00  -2.5 6.68e+00    -  1.00e+00 1.00e+00f  1
  91  6.6926290e+00 0.00e+00 2.37e+01  -2.5 2.24e+01    -  1.00e+00 5.00e-01f  2
  92  6.6920532e+00 0.00e+00 6.99e+00  -2.5 7.21e+00    -  1.00e+00 1.00e+00f  1
  93  6.6917366e+00 0.00e+00 2.17e+01  -2.5 2.17e+01    -  1.00e+00 5.00e-01f  2
  94  6.6912315e+00 0.00e+00 9.00e+00  -2.5 8.28e+00    -  1.00e+00 1.00e+00f  1
  95  6.6909337e+00 0.00e+00 1.68e+01  -2.5 1.86e+01    -  1.00e+00 5.00e-01f  2
  96  6.6905081e+00 0.00e+00 1.59e+01  -2.5 1.10e+01    -  1.00e+00 1.00e+00f  1
  97  6.6901136e+00 0.00e+00 1.83e+01  -2.5 1.19e+01    -  1.00e+00 1.00e+00f  1
  98  6.6897146e+00 0.00e+00 1.43e+01  -2.5 1.08e+01    -  1.00e+00 1.00e+00f  1
  99  6.6893860e+00 0.00e+00 2.39e+01  -2.5 1.39e+01    -  1.00e+00 1.00e+00f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
 100  6.6889879e+00 0.00e+00 8.99e+00  -2.5 9.11e+00    -  1.00e+00 1.00e+00f  1

Number of Iterations....: 100

                                   (scaled)                 (unscaled)
Objective...............:   4.1718859448364948e-01    6.6889879298493042e+00
Dual infeasibility......:   8.9941757858033142e+00    1.4420800104720774e+02
Constraint violation....:   0.0000000000000000e+00    0.0000000000000000e+00
Variable bound violation:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   0.0000000000000000e+00    0.0000000000000000e+00
Overall NLP error.......:   8.9941757858033142e+00    1.4420800104720774e+02


Number of objective function evaluations             = 300
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 seconds in IPOPT                               = 0.046

EXIT: Maximum Number of Iterations Exceeded.

┌ Info: * Solver : Ipopt
│ 
│ * Status
│   Result count       : 1
│   Termination status : ITERATION_LIMIT
│   Message from the solver:
│   "Maximum_Iterations_Exceeded"
│ 
│ * Candidate solution (result #1)
│   Primal status      : FEASIBLE_POINT
│   Dual status        : UNKNOWN_RESULT_STATUS
│   Objective value    : 6.68899e+00
│   Dual objective value : 0.00000e+00
│ 
│ * Work counters
│   Solve time (sec)   : 4.72462e-02
└   Barrier iterations : 100

(6.688987929849304, [384.54464659223163, 0.0005708194834985757, -377.29762681481145])

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
│   Result count       : 1
│   Termination status : LOCALLY_SOLVED
│   Message from the solver:
│   "Solve_Succeeded"
│ 
│ * Candidate solution (result #1)
│   Primal status      : FEASIBLE_POINT
│   Dual status        : FEASIBLE_POINT
│   Objective value    : 9.05568e-01
│   Dual objective value : 0.00000e+00
│ 
│ * Work counters
│   Solve time (sec)   : 6.61582e-01
└   Barrier iterations : 1745

(0.9055676718169392, [-1.9765000037247464, -0.6055995686432881, 8.481763938633279])
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é!

_, params = expfit(ts, ys, optimizer=Gurobi.Optimizer)

Ř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.