{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Applications of SOS Programming in Flowpipe Construction" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Flowpipe construction consists of under or over-approximating the sets of states reachable by dynamical systems. Recently a method has been developed for the class of polynomial ODEs with uncertain initial states (see [1], abbreviated `XFZ18under`). This method consists of reducing the Hamilton-Jacobi-Bellman equation to a hierarchy of semidefinite programs.\n", "\n", "In this notebook we consider the problem of approximating the flowpipe of a system of polynomial ODEs using `XFZ18under`. This is a Julia implementation that relies on the JuMP ecosystem (`JuMP`, `PolyJuMP`, `SumOfSquares`, `MathProgInterface`, `MathOptInterfaceMosek`) and the `JuliaAlgebra` ecosystem (`MultivariatePolynomials`, `DynamicPolynomials`). The implementation is evaluated on a set of standard benchmarks from formal verification and control engineering domains.\n", "\n", "---\n", "\n", "**References:**\n", "\n", "- [1] Xue, B., Fränzle, M., & Zhan, N. (2018, April). [Under-Approximating Reach Sets for Polynomial Continuous Systems. In Proceedings of the 21st International Conference on Hybrid Systems: Computation and Control (part of CPS Week) (pp. 51-60). ACM.](https://dl.acm.org/citation.cfm?id=3178133)\n", "\n", "**Quick links to documentation:**\n", "\n", "- www.juliaopt.org/SumOfSquares.jl/latest/\n", "- https://sums-of-squares.github.io/sos/index.html" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Van-der-Pol system\n", "\n", "Consider the following van-der-Pol system:\n", "\n", "$$\n", "\\dot{x}_1 = x_2 \\\\\n", "\\dot{x}_2 = -0.2x_1 + x_2 - 0.2x_1^2 x_2\n", "$$" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "differentiate (generic function with 18 methods)" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "using MultivariatePolynomials,\n", " JuMP,\n", " PolyJuMP,\n", " SumOfSquares,\n", " DynamicPolynomials,\n", " MathOptInterfaceMosek,\n", " MathematicalSystems,\n", " MosekTools\n", " \n", "const ∂ = differentiate" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/latex": [ "$$ ϵ $$" ], "text/plain": [ "ϵ" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# symbolic variables\n", "@polyvar x₁ x₂ t\n", "\n", "# time duration (scaled, see dynamics below)\n", "T = 1.0 \n", "\n", "# dynamics\n", "f = 2 * [x₂, -0.2*x₁ + x₂ - 0.2*x₁^2*x₂] \n", "\n", "# set of initial states X₀ = {x: V₀(x) <= 0}\n", "V₀ = x₁^2 + x₂^2 - 0.25\n", "\n", "# constraints Y = {x: g(x) >= 0} compact search space Y x [0, T]\n", "g = 25 - x₁^2 - x₂^2\n", "\n", "# degree of the relaxation\n", "k = 12\n", "\n", "# monomial vector up to order k, 0 <= sum_i alpha_i <= k, if alpha_i is the exponent of x_i\n", "X = monomials([x₁, x₂], 0:k)\n", "XT = monomials([x₁, x₂, t], 0:k)\n", "\n", "# create a SOS JuMP model to solve with Mosek\n", "model = SOSModel(with_optimizer(Mosek.Optimizer))\n", "\n", "# add unknown Φ to the model\n", "@variable(model, Φ, Poly(XT))\n", "\n", "# jacobian\n", "∂t = α -> ∂(α, t)\n", "∂xf = α -> ∂(α, x₁) * f[1] + ∂(α, x₂) * f[2] \n", "LΦ = ∂t(Φ) + ∂xf(Φ)\n", "\n", "# Φ(x, t) at time 0\n", "Φ₀ = subs(Φ, t => 0.)\n", "\n", "# scalar variable\n", "@variable(model, ϵ)\n", "\n", "dom1 = @set t*(T-t) >= 0 && g >= 0\n", "dom2 = @set g >= 0\n", "@constraint(model, ϵ >= 0.)\n", "@constraint(model, LΦ ∈ SOSCone(), domain = dom1)\n", "@constraint(model, ϵ - LΦ ∈ SOSCone(), domain = dom1)\n", "@constraint(model, Φ₀ - V₀ ∈ SOSCone(), domain = dom2)\n", "@constraint(model, ϵ + V₀ - Φ₀ ∈ SOSCone(), domain = dom2)\n", "\n", "@objective(model, Min, ϵ)" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Problem\n", " Name : \n", " Objective sense : min \n", " Type : CONIC (conic optimization problem)\n", " Constraints : 1543 \n", " Cones : 0 \n", " Scalar variables : 456 \n", " Matrix variables : 10 \n", " Integer variables : 0 \n", "\n", "Optimizer started.\n", "Presolve started.\n", "Linear dependency checker started.\n", "Linear dependency checker terminated.\n", "Eliminator started.\n", "Freed constraints in eliminator : 0\n", "Eliminator terminated.\n", "Eliminator - tries : 1 time : 0.00 \n", "Lin. dep. - tries : 1 time : 0.00 \n", "Lin. dep. - number : 0 \n", "Presolve terminated. Time: 0.00 \n", "Problem\n", " Name : \n", " Objective sense : min \n", " Type : CONIC (conic optimization problem)\n", " Constraints : 1543 \n", " Cones : 0 \n", " Scalar variables : 456 \n", " Matrix variables : 10 \n", " Integer variables : 0 \n", "\n", "Optimizer - threads : 8 \n", "Optimizer - solved problem : the primal \n", "Optimizer - Constraints : 1542\n", "Optimizer - Cones : 1\n", "Optimizer - Scalar variables : 457 conic : 456 \n", "Optimizer - Semi-definite variables: 10 scalarized : 30074 \n", "Factor - setup time : 0.16 dense det. time : 0.00 \n", "Factor - ML order time : 0.06 GP order time : 0.00 \n", "Factor - nonzeros before factor : 4.80e+05 after factor : 7.71e+05 \n", "Factor - dense dim. : 2 flops : 9.25e+08 \n", "ITE PFEAS DFEAS GFEAS PRSTATUS POBJ DOBJ MU TIME \n", "0 2.0e+00 8.0e+00 1.0e+00 0.00e+00 0.000000000e+00 0.000000000e+00 1.0e+00 0.18 \n", "1 9.1e-01 3.6e+00 6.2e-01 -7.85e-06 3.887272364e-02 -7.277377899e-04 4.6e-01 0.41 \n", "2 2.7e-01 1.1e+00 4.6e-01 9.08e-01 7.784730830e-01 4.303286309e-03 1.3e-01 0.65 \n", "3 6.2e-02 2.5e-01 4.5e-01 1.70e+00 2.542790179e-01 2.166139301e-02 3.1e-02 0.87 \n", "4 1.7e-02 6.7e-02 2.7e-01 1.87e+00 1.039863645e-01 6.446008794e-02 8.4e-03 1.10 \n", "5 9.8e-03 3.9e-02 2.1e-01 1.38e+00 1.294203141e-01 1.090778203e-01 4.9e-03 1.33 \n", "6 5.5e-03 2.2e-02 1.6e-01 1.15e+00 5.083213108e-02 4.002533190e-02 2.8e-03 1.55 \n", "7 1.9e-03 7.7e-03 1.0e-01 1.25e+00 2.537306998e-02 2.192283756e-02 9.7e-04 1.78 \n", "8 7.9e-04 3.2e-03 6.5e-02 1.18e+00 1.234485829e-02 1.103229702e-02 4.0e-04 2.00 \n", "9 1.7e-04 6.8e-04 3.2e-02 1.12e+00 3.993774249e-03 3.722721925e-03 8.5e-05 2.23 \n", "10 3.9e-05 1.6e-04 1.5e-02 1.04e+00 1.425668588e-03 1.364297146e-03 2.0e-05 2.45 \n", "11 1.1e-05 4.5e-05 7.1e-03 8.67e-01 8.002658418e-04 7.817030218e-04 5.7e-06 2.67 \n", "12 4.4e-06 1.8e-05 3.5e-03 6.29e-01 7.247429610e-04 7.167569170e-04 2.2e-06 2.89 \n", "13 1.1e-06 4.6e-06 9.2e-04 2.91e-01 1.130477537e-03 1.129671794e-03 5.7e-07 3.12 \n", "14 5.1e-07 2.0e-06 3.4e-04 -2.08e-01 1.970391763e-03 1.973234438e-03 2.5e-07 3.37 \n", "15 1.8e-07 6.0e-07 7.3e-05 -3.88e-01 4.513283862e-03 4.525581099e-03 7.3e-08 3.60 \n", "16 1.0e-07 1.7e-07 1.5e-05 -4.09e-01 9.367423708e-03 9.394061798e-03 2.0e-08 3.84 \n", "17 4.0e-08 6.3e-08 4.5e-06 -4.44e-01 1.682736656e-02 1.687283866e-02 7.6e-09 4.18 \n", "18 3.1e-08 3.2e-08 2.2e-06 -2.46e-01 2.305835380e-02 2.310891793e-02 3.8e-09 4.61 \n", "19 3.1e-08 3.2e-08 2.2e-06 -3.04e-01 2.305835380e-02 2.310891793e-02 3.8e-09 4.98 \n", "Optimizer terminated. Time: 5.25 \n", "\n", "Relaxation order : k = 12\n", "STALL\n", "JuMP.termination_status(model) = SLOW_PROGRESS\n", "JuMP.primal_status(model) = FEASIBLE_POINT\n", "JuMP.dual_status(model) = FEASIBLE_POINT\n", "STALL\n", "JuMP.objective_bound(model) = 0.0\n", "STALL\n", "JuMP.objective_value(model) = 0.02305835379792861\n" ] } ], "source": [ "optimize!(model)\n", "\n", "println(\"Relaxation order : k = $k\")\n", "println(\"JuMP.termination_status(model) = \", JuMP.termination_status(model))\n", "println(\"JuMP.primal_status(model) = \", JuMP.primal_status(model))\n", "println(\"JuMP.dual_status(model) = \", JuMP.dual_status(model))\n", "println(\"JuMP.objective_bound(model) = \", JuMP.objective_bound(model))\n", "println(\"JuMP.objective_value(model) = \", JuMP.objective_value(model))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Comparison with a YALMIP implementation in MATLAB\n", "\n", "| Package | k |Constraints|Scalar variables|Matrix variables|Time(s)|\n", "|-------------|------|-----------|----------------|----------------|-------|\n", "|JuMP v0.18.5 | 2 | 245 | 175 | 8 | < 1 |\n", "|**JuMP v0.19.0** | 2 | 83 | 13 | 8 | < 1 |\n", "|YALMIP | 2 | 152 | 63 | 10| < 1 |\n", "||\n", "|JuMP v0.18.5 | 3 | 893 | 715 | 10| ~ 1 |\n", "|**JuMP v0.19.0** | 3 | 199 | 21 | 10 | < 1 |\n", "|YALMIP | 3 | 254 | 121 | 10| ~ 1 |\n", "||\n", "|JuMP v0.18.5 | 4 | 893 | 730 | 10| ~1 |\n", "|**JuMP v0.19.0** | 4 | 199 | 36 | 10 | < 1 |\n", "|YALMIP | 4 | 394 | 206 | 10| 1.18 |\n", "||\n", "|JuMP v0.18.5 | 5 | 2639 | 2309 | 10 | 1.17 |\n", "|**JuMP v0.19.0** | 5 | 387 | 57 | 10 | < 1 |\n", "|YALMIP | 5 | 578 | 323 | 10 | 0.11 |\n", "||\n", "|JuMP v0.18.5 | 6 | 2639 | 2337 | 10| 0.90 |\n", "|**JuMP v0.19.0** | 6 | 387 | 85 | 10 | < 1 |\n", "|YALMIP | 6 | 812 | 477 | 10 | 1.10 |\n", "||\n", "|JuMP v0.18.5 | 7 | 6725 | 6183 | 10| 5.62 |\n", "|**JuMP v0.19.0** | 7 | 663 | 121 | 10 | < 1 |\n", "|YALMIP | 7 | 1102 | 673 | 10 | 1.52 |\n", "||\n", "|JuMP v0.18.5 | 8 | 6725 | 6228 | 10 | 6.06 |\n", "|**JuMP v0.19.0** | 8 | 663 | 166 | 10 | < 1 |\n", "|YALMIP | 8 | 1454 | 916 | 10 | 1.10 |\n", "||\n", "|JuMP v0.18.5 | 9 | 15269 | 14447 | 10 | 41.2 |\n", "|**JuMP v0.19.0** | 9 | 1043 | 221 | 10 | 1.70 |\n", "|YALMIP | 9 | 1874 | 1211 | 10 | 1.58 |\n", "||\n", "|JuMP v0.18.5 | 10 | 15269 | 14513 | 10 | 38.1 |\n", "|**JuMP v0.19.0** | 10 | 1043 | 287 | 10 | 1.67 |\n", "|YALMIP | 10 | 2368 | 1563 | 10 | 2.73 |\n", "||\n", "|JuMP v0.18.5 | 11 | 31617 | 30439 | 10 | 244 |\n", "|**JuMP v0.19.0** | 11 | 1543 | 365 | 10 | 4.88 |\n", "|YALMIP | 11 | 2942 | 1977 | 10 | 2.30 |\n", "||\n", "|JuMP v0.18.5 | 12 | 31617 | 30530 | 10 | 231 |\n", "|**JuMP v0.19.0** | 12 | 1543 | 456 | 10 | 5.02 |\n", "|YALMIP | 12 | 3602 | 2458 | 10 | 6.57 |\n", "\n", "Related: https://github.com/JuliaOpt/MosekTools.jl/issues/11." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "STALL\n" ] }, { "data": { "text/plain": [ "5.254102945327759" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "MOI.get(model, MOI.SolveTime())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Extracting the under and over approximations" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n" ] }, { "data": { "text/latex": [ "$$ -3.254406926195832e-6x₁^{12} + 3.834217164793429e-6x₁^{11}x₂ - 6.3384127954307375e-6x₁^{10}x₂^{2} + 1.7821158886139766e-6x₁^{9}x₂^{3} - 3.874093809241813e-7x₁^{8}x₂^{4} + 1.5483592082483843e-6x₁^{7}x₂^{5} - 1.7690055890869638e-6x₁^{6}x₂^{6} - 9.798790301107213e-7x₁^{5}x₂^{7} + 5.210797975903119e-7x₁^{4}x₂^{8} + 1.3500407479164792e-6x₁^{3}x₂^{9} + 4.4924259835334384e-7x₁^{2}x₂^{10} - 7.166251275333759e-8x₁x₂^{11} - 1.5032013847766846e-7x₂^{12} + 1.4831621204102143e-20x₁^{11} + 1.3691290786000017e-20x₁^{10}x₂ - 5.0093500950317737e-20x₁^{9}x₂^{2} + 4.3677787535910767e-20x₁^{8}x₂^{3} - 3.173767662076143e-20x₁^{7}x₂^{4} + 1.9724120886883362e-20x₁^{6}x₂^{5} + 3.932019772930437e-20x₁^{5}x₂^{6} - 3.366118784605938e-20x₁^{4}x₂^{7} + 7.335098065117538e-21x₁^{3}x₂^{8} + 1.8928666158525362e-21x₁^{2}x₂^{9} - 7.972721076141193e-21x₁x₂^{10} + 3.58418280872306e-21x₂^{11} + 0.00015991576845668712x₁^{10} - 0.0004657855037397151x₁^{9}x₂ + 0.000589085715209256x₁^{8}x₂^{2} - 0.0002325700371386792x₁^{7}x₂^{3} - 3.930451038774872e-5x₁^{6}x₂^{4} + 3.2322704287292744e-6x₁^{5}x₂^{5} + 1.3335992567483893e-6x₁^{4}x₂^{6} + 2.7518735707131795e-5x₁^{3}x₂^{7} - 1.323325367954451e-5x₁^{2}x₂^{8} - 1.2213506734162326e-5x₁x₂^{9} + 8.131677925472175e-6x₂^{10} - 2.873017741028833e-19x₁^{9} + 2.7918245968053696e-19x₁^{8}x₂ + 6.026007816753223e-19x₁^{7}x₂^{2} - 1.3364227498369648e-18x₁^{6}x₂^{3} + 1.01910474173647e-18x₁^{5}x₂^{4} - 1.8971371512328473e-19x₁^{4}x₂^{5} - 2.666149622126072e-19x₁^{3}x₂^{6} + 1.6020227232017553e-19x₁^{2}x₂^{7} - 2.0808839904916004e-20x₁x₂^{8} - 9.964099126024553e-21x₂^{9} - 0.0021077517585878706x₁^{8} + 0.007746530135541982x₁^{7}x₂ - 0.012729857429745972x₁^{6}x₂^{2} + 0.00959724703839421x₁^{5}x₂^{3} - 0.0026064147324748443x₁^{4}x₂^{4} - 0.0007538551940640615x₁^{3}x₂^{5} + 0.0008414955734000221x₁^{2}x₂^{6} - 0.0001935009339460489x₁x₂^{7} - 9.076112326494181e-6x₂^{8} + 3.0619458230114297e-18x₁^{7} - 6.825717234844004e-18x₁^{6}x₂ + 4.574693106727993e-18x₁^{5}x₂^{2} + 1.9761744010906326e-18x₁^{4}x₂^{3} - 5.1840687808885785e-18x₁^{3}x₂^{4} + 3.0469841339864345e-18x₁^{2}x₂^{5} - 6.108206321651024e-19x₁x₂^{6} + 7.73513798330052e-20x₂^{7} + 0.00766766678443157x₁^{6} - 0.02630247457713939x₁^{5}x₂ + 0.06439102517619352x₁^{4}x₂^{2} - 0.07522056249315182x₁^{3}x₂^{3} + 0.04442871000220708x₁^{2}x₂^{4} - 0.013324227892671911x₁x₂^{5} + 0.0018504331355650484x₂^{6} - 1.4946435441993517e-17x₁^{5} + 4.181493492797462e-17x₁^{4}x₂ - 5.630490775149093e-17x₁^{3}x₂^{2} + 4.678641430923887e-17x₁^{2}x₂^{3} - 2.2369497341344506e-17x₁x₂^{4} + 4.316626613003548e-18x₂^{5} - 0.026893175089150175x₁^{4} - 0.10307175783638588x₁^{3}x₂ + 0.23685478126822412x₁^{2}x₂^{2} - 0.14335941014008569x₁x₂^{3} + 0.033100700571824924x₂^{4} + 3.4341634754411844e-17x₁^{3} - 8.304800642136564e-17x₁^{2}x₂ + 7.510428581697817e-17x₁x₂^{2} - 2.3541145491650444e-17x₂^{3} + 0.6248870930138117x₁^{2} - 1.1487733659298884x₁x₂ + 0.5596157526120149x₂^{2} - 2.5712082025457398e-17x₁ + 2.3986229208007783e-17x₂ - 0.22309361984240117 $$" ], "text/plain": [ "-3.254406926195832e-6x₁¹² + 3.834217164793429e-6x₁¹¹x₂ - 6.3384127954307375e-6x₁¹⁰x₂² + 1.7821158886139766e-6x₁⁹x₂³ - 3.874093809241813e-7x₁⁸x₂⁴ + 1.5483592082483843e-6x₁⁷x₂⁵ - 1.7690055890869638e-6x₁⁶x₂⁶ - 9.798790301107213e-7x₁⁵x₂⁷ + 5.210797975903119e-7x₁⁴x₂⁸ + 1.3500407479164792e-6x₁³x₂⁹ + 4.4924259835334384e-7x₁²x₂¹⁰ - 7.166251275333759e-8x₁x₂¹¹ - 1.5032013847766846e-7x₂¹² + 1.4831621204102143e-20x₁¹¹ + 1.3691290786000017e-20x₁¹⁰x₂ - 5.0093500950317737e-20x₁⁹x₂² + 4.3677787535910767e-20x₁⁸x₂³ - 3.173767662076143e-20x₁⁷x₂⁴ + 1.9724120886883362e-20x₁⁶x₂⁵ + 3.932019772930437e-20x₁⁵x₂⁶ - 3.366118784605938e-20x₁⁴x₂⁷ + 7.335098065117538e-21x₁³x₂⁸ + 1.8928666158525362e-21x₁²x₂⁹ - 7.972721076141193e-21x₁x₂¹⁰ + 3.58418280872306e-21x₂¹¹ + 0.00015991576845668712x₁¹⁰ - 0.0004657855037397151x₁⁹x₂ + 0.000589085715209256x₁⁸x₂² - 0.0002325700371386792x₁⁷x₂³ - 3.930451038774872e-5x₁⁶x₂⁴ + 3.2322704287292744e-6x₁⁵x₂⁵ + 1.3335992567483893e-6x₁⁴x₂⁶ + 2.7518735707131795e-5x₁³x₂⁷ - 1.323325367954451e-5x₁²x₂⁸ - 1.2213506734162326e-5x₁x₂⁹ + 8.131677925472175e-6x₂¹⁰ - 2.873017741028833e-19x₁⁹ + 2.7918245968053696e-19x₁⁸x₂ + 6.026007816753223e-19x₁⁷x₂² - 1.3364227498369648e-18x₁⁶x₂³ + 1.01910474173647e-18x₁⁵x₂⁴ - 1.8971371512328473e-19x₁⁴x₂⁵ - 2.666149622126072e-19x₁³x₂⁶ + 1.6020227232017553e-19x₁²x₂⁷ - 2.0808839904916004e-20x₁x₂⁸ - 9.964099126024553e-21x₂⁹ - 0.0021077517585878706x₁⁸ + 0.007746530135541982x₁⁷x₂ - 0.012729857429745972x₁⁶x₂² + 0.00959724703839421x₁⁵x₂³ - 0.0026064147324748443x₁⁴x₂⁴ - 0.0007538551940640615x₁³x₂⁵ + 0.0008414955734000221x₁²x₂⁶ - 0.0001935009339460489x₁x₂⁷ - 9.076112326494181e-6x₂⁸ + 3.0619458230114297e-18x₁⁷ - 6.825717234844004e-18x₁⁶x₂ + 4.574693106727993e-18x₁⁵x₂² + 1.9761744010906326e-18x₁⁴x₂³ - 5.1840687808885785e-18x₁³x₂⁴ + 3.0469841339864345e-18x₁²x₂⁵ - 6.108206321651024e-19x₁x₂⁶ + 7.73513798330052e-20x₂⁷ + 0.00766766678443157x₁⁶ - 0.02630247457713939x₁⁵x₂ + 0.06439102517619352x₁⁴x₂² - 0.07522056249315182x₁³x₂³ + 0.04442871000220708x₁²x₂⁴ - 0.013324227892671911x₁x₂⁵ + 0.0018504331355650484x₂⁶ - 1.4946435441993517e-17x₁⁵ + 4.181493492797462e-17x₁⁴x₂ - 5.630490775149093e-17x₁³x₂² + 4.678641430923887e-17x₁²x₂³ - 2.2369497341344506e-17x₁x₂⁴ + 4.316626613003548e-18x₂⁵ - 0.026893175089150175x₁⁴ - 0.10307175783638588x₁³x₂ + 0.23685478126822412x₁²x₂² - 0.14335941014008569x₁x₂³ + 0.033100700571824924x₂⁴ + 3.4341634754411844e-17x₁³ - 8.304800642136564e-17x₁²x₂ + 7.510428581697817e-17x₁x₂² - 2.3541145491650444e-17x₂³ + 0.6248870930138117x₁² - 1.1487733659298884x₁x₂ + 0.5596157526120149x₂² - 2.5712082025457398e-17x₁ + 2.3986229208007783e-17x₂ - 0.22309361984240117" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Recovering the solution:\n", "ϵopt = JuMP.objective_value(model)\n", "\n", "# Punder <= 0\n", "Punder = subs(JuMP.value(model[:Φ]), t => T)" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n", "STALL\n" ] }, { "data": { "text/latex": [ "$$ -3.254406926195832e-6x₁^{12} + 3.834217164793429e-6x₁^{11}x₂ - 6.3384127954307375e-6x₁^{10}x₂^{2} + 1.7821158886139766e-6x₁^{9}x₂^{3} - 3.874093809241813e-7x₁^{8}x₂^{4} + 1.5483592082483843e-6x₁^{7}x₂^{5} - 1.7690055890869638e-6x₁^{6}x₂^{6} - 9.798790301107213e-7x₁^{5}x₂^{7} + 5.210797975903119e-7x₁^{4}x₂^{8} + 1.3500407479164792e-6x₁^{3}x₂^{9} + 4.4924259835334384e-7x₁^{2}x₂^{10} - 7.166251275333759e-8x₁x₂^{11} - 1.5032013847766846e-7x₂^{12} + 1.4831621204102143e-20x₁^{11} + 1.3691290786000017e-20x₁^{10}x₂ - 5.0093500950317737e-20x₁^{9}x₂^{2} + 4.3677787535910767e-20x₁^{8}x₂^{3} - 3.173767662076143e-20x₁^{7}x₂^{4} + 1.9724120886883362e-20x₁^{6}x₂^{5} + 3.932019772930437e-20x₁^{5}x₂^{6} - 3.366118784605938e-20x₁^{4}x₂^{7} + 7.335098065117538e-21x₁^{3}x₂^{8} + 1.8928666158525362e-21x₁^{2}x₂^{9} - 7.972721076141193e-21x₁x₂^{10} + 3.58418280872306e-21x₂^{11} + 0.00015991576845668712x₁^{10} - 0.0004657855037397151x₁^{9}x₂ + 0.000589085715209256x₁^{8}x₂^{2} - 0.0002325700371386792x₁^{7}x₂^{3} - 3.930451038774872e-5x₁^{6}x₂^{4} + 3.2322704287292744e-6x₁^{5}x₂^{5} + 1.3335992567483893e-6x₁^{4}x₂^{6} + 2.7518735707131795e-5x₁^{3}x₂^{7} - 1.323325367954451e-5x₁^{2}x₂^{8} - 1.2213506734162326e-5x₁x₂^{9} + 8.131677925472175e-6x₂^{10} - 2.873017741028833e-19x₁^{9} + 2.7918245968053696e-19x₁^{8}x₂ + 6.026007816753223e-19x₁^{7}x₂^{2} - 1.3364227498369648e-18x₁^{6}x₂^{3} + 1.01910474173647e-18x₁^{5}x₂^{4} - 1.8971371512328473e-19x₁^{4}x₂^{5} - 2.666149622126072e-19x₁^{3}x₂^{6} + 1.6020227232017553e-19x₁^{2}x₂^{7} - 2.0808839904916004e-20x₁x₂^{8} - 9.964099126024553e-21x₂^{9} - 0.0021077517585878706x₁^{8} + 0.007746530135541982x₁^{7}x₂ - 0.012729857429745972x₁^{6}x₂^{2} + 0.00959724703839421x₁^{5}x₂^{3} - 0.0026064147324748443x₁^{4}x₂^{4} - 0.0007538551940640615x₁^{3}x₂^{5} + 0.0008414955734000221x₁^{2}x₂^{6} - 0.0001935009339460489x₁x₂^{7} - 9.076112326494181e-6x₂^{8} + 3.0619458230114297e-18x₁^{7} - 6.825717234844004e-18x₁^{6}x₂ + 4.574693106727993e-18x₁^{5}x₂^{2} + 1.9761744010906326e-18x₁^{4}x₂^{3} - 5.1840687808885785e-18x₁^{3}x₂^{4} + 3.0469841339864345e-18x₁^{2}x₂^{5} - 6.108206321651024e-19x₁x₂^{6} + 7.73513798330052e-20x₂^{7} + 0.00766766678443157x₁^{6} - 0.02630247457713939x₁^{5}x₂ + 0.06439102517619352x₁^{4}x₂^{2} - 0.07522056249315182x₁^{3}x₂^{3} + 0.04442871000220708x₁^{2}x₂^{4} - 0.013324227892671911x₁x₂^{5} + 0.0018504331355650484x₂^{6} - 1.4946435441993517e-17x₁^{5} + 4.181493492797462e-17x₁^{4}x₂ - 5.630490775149093e-17x₁^{3}x₂^{2} + 4.678641430923887e-17x₁^{2}x₂^{3} - 2.2369497341344506e-17x₁x₂^{4} + 4.316626613003548e-18x₂^{5} - 0.026893175089150175x₁^{4} - 0.10307175783638588x₁^{3}x₂ + 0.23685478126822412x₁^{2}x₂^{2} - 0.14335941014008569x₁x₂^{3} + 0.033100700571824924x₂^{4} + 3.4341634754411844e-17x₁^{3} - 8.304800642136564e-17x₁^{2}x₂ + 7.510428581697817e-17x₁x₂^{2} - 2.3541145491650444e-17x₂^{3} + 0.6248870930138117x₁^{2} - 1.1487733659298884x₁x₂ + 0.5596157526120149x₂^{2} - 2.5712082025457398e-17x₁ + 2.3986229208007783e-17x₂ - 0.2692103274382584 $$" ], "text/plain": [ "-3.254406926195832e-6x₁¹² + 3.834217164793429e-6x₁¹¹x₂ - 6.3384127954307375e-6x₁¹⁰x₂² + 1.7821158886139766e-6x₁⁹x₂³ - 3.874093809241813e-7x₁⁸x₂⁴ + 1.5483592082483843e-6x₁⁷x₂⁵ - 1.7690055890869638e-6x₁⁶x₂⁶ - 9.798790301107213e-7x₁⁵x₂⁷ + 5.210797975903119e-7x₁⁴x₂⁸ + 1.3500407479164792e-6x₁³x₂⁹ + 4.4924259835334384e-7x₁²x₂¹⁰ - 7.166251275333759e-8x₁x₂¹¹ - 1.5032013847766846e-7x₂¹² + 1.4831621204102143e-20x₁¹¹ + 1.3691290786000017e-20x₁¹⁰x₂ - 5.0093500950317737e-20x₁⁹x₂² + 4.3677787535910767e-20x₁⁸x₂³ - 3.173767662076143e-20x₁⁷x₂⁴ + 1.9724120886883362e-20x₁⁶x₂⁵ + 3.932019772930437e-20x₁⁵x₂⁶ - 3.366118784605938e-20x₁⁴x₂⁷ + 7.335098065117538e-21x₁³x₂⁸ + 1.8928666158525362e-21x₁²x₂⁹ - 7.972721076141193e-21x₁x₂¹⁰ + 3.58418280872306e-21x₂¹¹ + 0.00015991576845668712x₁¹⁰ - 0.0004657855037397151x₁⁹x₂ + 0.000589085715209256x₁⁸x₂² - 0.0002325700371386792x₁⁷x₂³ - 3.930451038774872e-5x₁⁶x₂⁴ + 3.2322704287292744e-6x₁⁵x₂⁵ + 1.3335992567483893e-6x₁⁴x₂⁶ + 2.7518735707131795e-5x₁³x₂⁷ - 1.323325367954451e-5x₁²x₂⁸ - 1.2213506734162326e-5x₁x₂⁹ + 8.131677925472175e-6x₂¹⁰ - 2.873017741028833e-19x₁⁹ + 2.7918245968053696e-19x₁⁸x₂ + 6.026007816753223e-19x₁⁷x₂² - 1.3364227498369648e-18x₁⁶x₂³ + 1.01910474173647e-18x₁⁵x₂⁴ - 1.8971371512328473e-19x₁⁴x₂⁵ - 2.666149622126072e-19x₁³x₂⁶ + 1.6020227232017553e-19x₁²x₂⁷ - 2.0808839904916004e-20x₁x₂⁸ - 9.964099126024553e-21x₂⁹ - 0.0021077517585878706x₁⁸ + 0.007746530135541982x₁⁷x₂ - 0.012729857429745972x₁⁶x₂² + 0.00959724703839421x₁⁵x₂³ - 0.0026064147324748443x₁⁴x₂⁴ - 0.0007538551940640615x₁³x₂⁵ + 0.0008414955734000221x₁²x₂⁶ - 0.0001935009339460489x₁x₂⁷ - 9.076112326494181e-6x₂⁸ + 3.0619458230114297e-18x₁⁷ - 6.825717234844004e-18x₁⁶x₂ + 4.574693106727993e-18x₁⁵x₂² + 1.9761744010906326e-18x₁⁴x₂³ - 5.1840687808885785e-18x₁³x₂⁴ + 3.0469841339864345e-18x₁²x₂⁵ - 6.108206321651024e-19x₁x₂⁶ + 7.73513798330052e-20x₂⁷ + 0.00766766678443157x₁⁶ - 0.02630247457713939x₁⁵x₂ + 0.06439102517619352x₁⁴x₂² - 0.07522056249315182x₁³x₂³ + 0.04442871000220708x₁²x₂⁴ - 0.013324227892671911x₁x₂⁵ + 0.0018504331355650484x₂⁶ - 1.4946435441993517e-17x₁⁵ + 4.181493492797462e-17x₁⁴x₂ - 5.630490775149093e-17x₁³x₂² + 4.678641430923887e-17x₁²x₂³ - 2.2369497341344506e-17x₁x₂⁴ + 4.316626613003548e-18x₂⁵ - 0.026893175089150175x₁⁴ - 0.10307175783638588x₁³x₂ + 0.23685478126822412x₁²x₂² - 0.14335941014008569x₁x₂³ + 0.033100700571824924x₂⁴ + 3.4341634754411844e-17x₁³ - 8.304800642136564e-17x₁²x₂ + 7.510428581697817e-17x₁x₂² - 2.3541145491650444e-17x₂³ + 0.6248870930138117x₁² - 1.1487733659298884x₁x₂ + 0.5596157526120149x₂² - 2.5712082025457398e-17x₁ + 2.3986229208007783e-17x₂ - 0.2692103274382584" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Pover <= 0\n", "Pover = subs(JuMP.value(model[:Φ]), t => T) - ϵopt * (T+1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Visualization using IntervalConstraintProgramming\n", "\n", "See notebook `Constraints.ipynb`." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# no necesito convertir, si utilizo el IntervalContractors#fix_zero_power\n", "using ModelingToolkit\n", "\n", "# TODO: temporary definition here, needs convert(Variable, z^2 + t^2)\n", "vars = ModelingToolkit.@variables x1 x2\n", "\n", "function tomtexpr(poly::Polynomial)\n", " global vars\n", " expr = 0.0\n", " for (i, term) in enumerate(poly.x)\n", " if iszero(exponents(term))\n", " expr += poly.a[i]\n", " else\n", " expr += poly.a[i] * prod([vars[j]^value for (j, value) in enumerate(exponents(term)) if value != 0])\n", " end\n", " end\n", " return expr\n", "end" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "Punder_mt = tomtexpr(Punder)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "include(\"constraints.jl\")\n", "\n", "using Plots" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "S = Separator(-Inf..0.0, Punder_mt)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "p = pave(S, IntervalBox(-2..2, 2), 0.1)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "S = Separator(-Inf..0.0, Punder_mt)\n", "p = pave(S, IntervalBox(-3..3, 2), 0.1)\n", "\n", "plot(p.inner, lw=0, aspectratio=1, label=\"inner\")\n", "plot!(p.boundary, lw=0, label=\"boundary\", legend=:topleft)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "Punder(2.0, 0.5)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "S = Separator(-Inf..0.0, Punder_mt)\n", "p = pave(S, IntervalBox(-2..2, 2), 0.1)\n", "\n", "Pover_mt = tomtexpr(Pover)\n", "Sover = Separator(-Inf..0.0, Pover_mt)\n", "pover = pave(Sover, IntervalBox(-2..2, 2), 0.1)\n", "\n", "plot(p.inner, lw=0, aspectratio=1, label=\"under\", alpha=.4)\n", "plot!(pover.inner, lw=0, label=\"over\", legend=:topleft, alpha=.4)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Visualization" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "using LinearAlgebra, StaticPolynomials\n", "\n", "# filter out the smallest coefficients?\n", "filter_coeffs = false\n", "\n", "if filter_coeffs\n", " max_coeff = maximum(abs.(Punder.a))\n", "\n", " # plots are *very* sensitive to small values\n", " # make tol bigger than 0 to filter out some coefficients\n", " tol = 0.0 # max_coeff / 1e19\n", " kept_coeffs = map(c -> abs(c) > tol, Punder.a)\n", " Punder = dot(Punder.a[kept_coeffs], Punder.x[kept_coeffs])\n", " Pover = dot(Pover.a[kept_coeffs], Pover.x[kept_coeffs]);\n", "end\n", "\n", "# we now convert to a static polynomial for faster evaluation\n", "Punder_st = StaticPolynomials.Polynomial(Punder)\n", "Pover_st = StaticPolynomials.Polynomial(Pover);\n", "\n", "_Punder(x, y) = Punder_st([x, y])\n", "_Pover(x, y) = Pover_st([x, y])" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "@btime Punder(1.0, 1.0)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "@btime Punder_st([1.0, 1.0])" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "2e-6 / 96e-9" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "Punder_st" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Plots using ImplicitEquations" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "using ImplicitEquations, Plots\n", "gr()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "G = plot()\n", "plot!(G, _Punder ⩵ 0., xlims=(-4, 4), ylims=(-4, 4), color=\"red\")\n", "plot!(G, _Pover ⩵ 0., xlims=(-4, 4), ylims=(-4, 4), color=\"blue\")\n", "G" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "G = plot()\n", "Gu = plot(_Punder ≪ 0., xlims=(-2, 2), ylims=(-2, 2))\n", "Go = plot(_Pover ≪ 0., xlims=(-2, 2), ylims=(-2, 2))\n", "plot(Gu, Go)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "Punder(2.0, 1.5)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Plots using IntervalConstraintProgramming" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "using IntervalConstraintProgramming, ValidatedNumerics" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "@function Funder(x₁, x₂) = _Punder(x₁, x₂)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "Bx₁x₂ = IntervalBox(-10..10, -10..10)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "Cunder = IntervalConstraintProgramming.@constraint Funder(x₁, x₂) <= 0.0 # no" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "paving = pave(Cunder, Bx₁x₂)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Adding a constraint" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "text/latex": [ "$$ ϵ $$" ], "text/plain": [ "ϵ" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# symbolic variables\n", "@polyvar x₁ x₂ t\n", "\n", "# time duration (scaled, see dynamics below)\n", "T = 1.0 \n", "\n", "# dynamics\n", "f = 2 * [x₂, -0.2*x₁ + x₂ - 0.2*x₁^2*x₂] \n", "\n", "# set of initial states X₀ = {x: V₀(x) <= 0}\n", "V₀ = x₁^2 + x₂^2 - 0.25\n", "\n", "# constraints Y = {x: g(x) >= 0} compact search space Y x [0, T]\n", "g = 25 - x₁^2 - x₂^2\n", "\n", "# degree of the relaxation\n", "k = 12\n", "\n", "# monomial vector up to order k, 0 <= sum_i alpha_i <= k, if alpha_i is the exponent of x_i\n", "X = monomials([x₁, x₂], 0:k)\n", "XT = monomials([x₁, x₂, t], 0:k)\n", "\n", "# create a SOS JuMP model to solve with Mosek\n", "model = SOSModel(with_optimizer(Mosek.Optimizer))\n", "\n", "# add unknown Φ to the model\n", "@variable(model, Φ, Poly(XT))\n", "\n", "# jacobian\n", "∂t = α -> ∂(α, t)\n", "∂xf = α -> ∂(α, x₁) * f[1] + ∂(α, x₂) * f[2] \n", "LΦ = ∂t(Φ) + ∂xf(Φ)\n", "\n", "# Φ(x, t) at time 0\n", "Φ₀ = subs(Φ, t => 0.)\n", "\n", "# scalar variable\n", "@variable(model, ϵ)\n", "\n", "dom1 = @set t*(T-t) >= 0 && g >= 0\n", "dom2 = @set g >= 0\n", "@constraint(model, ϵ >= 0.)\n", "@constraint(model, LΦ ∈ SOSCone(), domain = dom1)\n", "@constraint(model, ϵ - LΦ ∈ SOSCone(), domain = dom1)\n", "@constraint(model, Φ₀ - V₀ ∈ SOSCone(), domain = dom2)\n", "@constraint(model, ϵ + V₀ - Φ₀ ∈ SOSCone(), domain = dom2)\n", "\n", "# add a safety constarint\n", "#@constraint(model, x₂ <= 100.)\n", "\n", "@objective(model, Min, ϵ)" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Problem\n", " Name : \n", " Objective sense : min \n", " Type : CONIC (conic optimization problem)\n", " Constraints : 1543 \n", " Cones : 0 \n", " Scalar variables : 456 \n", " Matrix variables : 10 \n", " Integer variables : 0 \n", "\n", "Optimizer started.\n", "Presolve started.\n", "Linear dependency checker started.\n", "Linear dependency checker terminated.\n", "Eliminator started.\n", "Freed constraints in eliminator : 0\n", "Eliminator terminated.\n", "Eliminator - tries : 1 time : 0.00 \n", "Lin. dep. - tries : 1 time : 0.00 \n", "Lin. dep. - number : 0 \n", "Presolve terminated. Time: 0.00 \n", "Problem\n", " Name : \n", " Objective sense : min \n", " Type : CONIC (conic optimization problem)\n", " Constraints : 1543 \n", " Cones : 0 \n", " Scalar variables : 456 \n", " Matrix variables : 10 \n", " Integer variables : 0 \n", "\n", "Optimizer - threads : 8 \n", "Optimizer - solved problem : the primal \n", "Optimizer - Constraints : 1542\n", "Optimizer - Cones : 1\n", "Optimizer - Scalar variables : 457 conic : 456 \n", "Optimizer - Semi-definite variables: 10 scalarized : 30074 \n", "Factor - setup time : 0.15 dense det. time : 0.00 \n", "Factor - ML order time : 0.06 GP order time : 0.00 \n", "Factor - nonzeros before factor : 4.80e+05 after factor : 7.71e+05 \n", "Factor - dense dim. : 2 flops : 9.25e+08 \n", "ITE PFEAS DFEAS GFEAS PRSTATUS POBJ DOBJ MU TIME \n", "0 2.0e+00 8.0e+00 1.0e+00 0.00e+00 0.000000000e+00 0.000000000e+00 1.0e+00 0.17 \n", "1 9.1e-01 3.6e+00 6.2e-01 -7.85e-06 3.887272364e-02 -7.277377899e-04 4.6e-01 0.39 \n", "2 2.7e-01 1.1e+00 4.6e-01 9.08e-01 7.784730830e-01 4.303286309e-03 1.3e-01 0.61 \n", "3 6.2e-02 2.5e-01 4.5e-01 1.70e+00 2.542790179e-01 2.166139301e-02 3.1e-02 0.81 \n", "4 1.7e-02 6.7e-02 2.7e-01 1.87e+00 1.039863645e-01 6.446008794e-02 8.4e-03 1.08 \n", "5 9.8e-03 3.9e-02 2.1e-01 1.38e+00 1.294203141e-01 1.090778203e-01 4.9e-03 1.30 \n", "6 5.5e-03 2.2e-02 1.6e-01 1.15e+00 5.083213108e-02 4.002533190e-02 2.8e-03 1.54 \n", "7 1.9e-03 7.7e-03 1.0e-01 1.25e+00 2.537306998e-02 2.192283756e-02 9.7e-04 1.78 \n", "8 7.9e-04 3.2e-03 6.5e-02 1.18e+00 1.234485829e-02 1.103229702e-02 4.0e-04 2.01 \n", "9 1.7e-04 6.8e-04 3.2e-02 1.12e+00 3.993774249e-03 3.722721925e-03 8.5e-05 2.22 \n", "10 3.9e-05 1.6e-04 1.5e-02 1.04e+00 1.425668588e-03 1.364297146e-03 2.0e-05 2.48 \n", "11 1.1e-05 4.5e-05 7.1e-03 8.67e-01 8.002658418e-04 7.817030218e-04 5.7e-06 2.69 \n", "12 4.4e-06 1.8e-05 3.5e-03 6.29e-01 7.247429610e-04 7.167569170e-04 2.2e-06 2.94 \n", "13 1.1e-06 4.6e-06 9.2e-04 2.91e-01 1.130477537e-03 1.129671794e-03 5.7e-07 3.15 \n", "14 5.1e-07 2.0e-06 3.4e-04 -2.08e-01 1.970391763e-03 1.973234438e-03 2.5e-07 3.38 \n", "15 1.8e-07 6.0e-07 7.3e-05 -3.88e-01 4.513283862e-03 4.525581099e-03 7.3e-08 3.58 \n", "16 1.0e-07 1.7e-07 1.5e-05 -4.09e-01 9.367423708e-03 9.394061798e-03 2.0e-08 3.79 \n", "17 4.0e-08 6.3e-08 4.5e-06 -4.44e-01 1.682736656e-02 1.687283866e-02 7.6e-09 4.08 \n", "18 3.1e-08 3.2e-08 2.2e-06 -2.46e-01 2.305835380e-02 2.310891793e-02 3.8e-09 4.48 \n", "19 3.1e-08 3.2e-08 2.2e-06 -3.04e-01 2.305835380e-02 2.310891793e-02 3.8e-09 4.83 \n", "Optimizer terminated. Time: 5.01 \n", "\n", "Relaxation order : k = 12\n", "STALL\n", "JuMP.termination_status(model) = SLOW_PROGRESS\n", "JuMP.primal_status(model) = FEASIBLE_POINT\n", "JuMP.dual_status(model) = FEASIBLE_POINT\n", "STALL\n", "JuMP.objective_bound(model) = 0.0\n", "STALL\n", "JuMP.objective_value(model) = 0.02305835379792861\n" ] } ], "source": [ "optimize!(model)\n", "\n", "println(\"Relaxation order : k = $k\")\n", "println(\"JuMP.termination_status(model) = \", JuMP.termination_status(model))\n", "println(\"JuMP.primal_status(model) = \", JuMP.primal_status(model))\n", "println(\"JuMP.dual_status(model) = \", JuMP.dual_status(model))\n", "println(\"JuMP.objective_bound(model) = \", JuMP.objective_bound(model))\n", "println(\"JuMP.objective_value(model) = \", JuMP.objective_value(model))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# symbolic variables\n", "@polyvar x₁ x₂ t\n", "\n", "# time duration (scaled, see dynamics below)\n", "T = 1.0 \n", "\n", "# dynamics\n", "f = 2 * [x₂, -0.2*x₁ + x₂ - 0.2*x₁^2*x₂] \n", "\n", "# set of initial states X₀ = {x: V₀(x) <= 0}\n", "V₀ = x₁^2 + x₂^2 - 0.25\n", "\n", "# constraints Y = {x: g(x) >= 0} compact search space Y x [0, T]\n", "g = 25 - x₁^2 - x₂^2\n", "\n", "# degree of the relaxation\n", "k = 12\n", "\n", "# monomial vector up to order k, 0 <= sum_i alpha_i <= k, if alpha_i is the exponent of x_i\n", "X = monomials([x₁, x₂], 0:k)\n", "XT = monomials([x₁, x₂, t], 0:k)\n", "\n", "# create a SOS JuMP model to solve with Mosek\n", "model_verif = SOSModel(with_optimizer(Mosek.Optimizer))\n", "\n", "# add unknown Φ to the model\n", "@variable(model, Φ, Poly(XT))\n", "\n", "# jacobian\n", "∂t = α -> ∂(α, t)\n", "∂xf = α -> ∂(α, x₁) * f[1] + ∂(α, x₂) * f[2] \n", "LΦ = ∂t(Φ) + ∂xf(Φ)\n", "\n", "# Φ(x, t) at time 0\n", "Φ₀ = subs(Φ, t => 0.)\n", "\n", "# scalar variable\n", "@variable(model, ϵ)\n", "\n", "dom1 = @set t*(T-t) >= 0 && g >= 0\n", "dom2 = @set g >= 0\n", "@constraint(model, ϵ >= 0.)\n", "@constraint(model, LΦ ∈ SOSCone(), domain = dom1)\n", "@constraint(model, ϵ - LΦ ∈ SOSCone(), domain = dom1)\n", "@constraint(model, Φ₀ - V₀ ∈ SOSCone(), domain = dom2)\n", "@constraint(model, ϵ + V₀ - Φ₀ ∈ SOSCone(), domain = dom2)\n", "\n", "# add a safety constarint\n", "#@constraint(model, x₂ <= 100.)\n", "\n", "@objective(model, Min, ϵ)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "Punder" ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [ { "ename": "ErrorException", "evalue": "PolyJuMP is just a JuMP extension for modelling Polynomial Optimization: it does not implement any reformulation. To use automatic sums of squares (SOS) reformulations, install the SumOfSquares Julia package and try `using SumOfSquares` and `setpolymodule!(SumOfSquares)` or use `SOSModel` instead of `Model`.", "output_type": "error", "traceback": [ "PolyJuMP is just a JuMP extension for modelling Polynomial Optimization: it does not implement any reformulation. To use automatic sums of squares (SOS) reformulations, install the SumOfSquares Julia package and try `using SumOfSquares` and `setpolymodule!(SumOfSquares)` or use `SOSModel` instead of `Model`.", "", "Stacktrace:", " [1] error(::String) at ./error.jl:33", " [2] get_default(::Nothing) at /Users/forets/.julia/packages/PolyJuMP/DffqD/src/default.jl:20", " [3] getdefault(::PolyJuMP.Data, ::Type{PolyJuMP.NonNegPoly}) at /Users/forets/.julia/packages/PolyJuMP/DffqD/src/default.jl:12", " [4] getdefault(::PolyJuMP.Data, ::PolyJuMP.NonNegPoly) at /Users/forets/.julia/packages/PolyJuMP/DffqD/src/default.jl:11", " [5] getdefault(::Model, ::PolyJuMP.NonNegPoly) at /Users/forets/.julia/packages/PolyJuMP/DffqD/src/default.jl:7", " [6] add_constraint(::Model, ::PolyJuMP.Constraint{getfield(JuMP, Symbol(\"#_error#55\")){Symbol},Polynomial{true,Float64},PolyJuMP.NonNegPoly,Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}}, ::String) at /Users/forets/.julia/packages/PolyJuMP/DffqD/src/constraint.jl:166", " [7] top-level scope at /Users/forets/.julia/packages/JuMP/jnmGG/src/macros.jl:621", " [8] top-level scope at In[29]:6" ] } ], "source": [ "model_verif = Model(with_optimizer(Mosek.Optimizer))\n", "\n", "@variable(model_verif, t)\n", "@variable(model_verif, x₁)\n", "@variable(model_verif, x₂)\n", "\n", "@constraint(model_verif, Punder <= 0)\n", "\n", "@objective(model_verif, Max, x₂)\n", "\n", "optimize!(model_verif)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Julia 1.1.0-rc2", "language": "julia", "name": "julia-1.1" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.1.0" } }, "nbformat": 4, "nbformat_minor": 2 }