### A Pluto.jl notebook ###
# v0.19.46
using Markdown
using InteractiveUtils
# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
quote
local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end
local el = $(esc(element))
global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
el
end
end
# ╔═╡ 9c61c9d5-0d3d-47b7-8e7d-f929cfec3311
begin
#
using LaTeXStrings
using JuMP
using Plots
using Ipopt
using GLPK
using Distributions
using Statistics
using OrderedCollections
using LinearAlgebra
using Test
using NLopt
using PlutoUI
end
# ╔═╡ aa1f0719-5e1d-4294-815e-349d1cebc004
html"""
"""
# ╔═╡ aae2e545-198c-42bb-9e97-6d6635bee02a
html""
# ╔═╡ 4468041e-4472-11ec-16b7-553f4def2b18
md"""
# Constraints
Recall our core optimization problem:
$$\min_{x\in\mathbb{R}^n} f(x) \text{ s.t. } x \in \mathcal{X}$$
* Up to now, the feasible set was $\mathcal{X} \in \mathbb{R}^n$.
* In **constrained problems** $\mathcal{X}$ is a subset thereof.
* We already encountered *box constraints*, e.g. $x \in [a,b]$.
* Sometimes the contrained solution coincides with the unconstrained one, sometimes it does not.
* There are *equality constraints* and *inequality constraints*.
"""
# ╔═╡ e66005ad-1f51-462a-9329-d77f7695cbfc
md"""
## Lagrange Multipliers
* Used to optimize a function subject to equality constraints.
$$\begin{aligned}
\min_x & f(x) \\
\text{subject to } & h(x) = 0
\end{aligned}$$
where both $f$ and $h$ have continuous partial derivatives.
* We look for contour lines of $f$ that are aligned to contours of $h(x) = 0$.
In other words, we want to find the best $x$ s.t. $h(x) = 0$ and we have
$$\nabla f(x) = \lambda \nabla h(x)$$
for some *Lagrange Mutliplier* $\lambda$
* Notice that we need the scalar $\lambda$ because the magnitudes of the gradients may be different.
* We therefore form the the **Lagrangian**:
$$\mathcal{L}(x,\lambda) = f(x) - \lambda h(x)$$
"""
# ╔═╡ 5bc7c8ce-3a4a-4505-a2fa-2b25f041437a
md"""
## Example
Suppose we have
$$\begin{align}
\min_x & -\exp\left( -\left( x_1 x_2 - \frac{3}{2} \right)^2 - \left(x_2 - \frac{3}{2}\right)^2 \right) \\
\text{subject to } & x_1 - x_2^2 = 0
\end{align}$$
We form the Lagrangiagn:
$$\mathcal{L}(x_1,x_2,\lambda) = -\exp\left( -\left( x_1 x_2 - \frac{3}{2} \right)^2 - \left(x_2 - \frac{3}{2}\right)^2 \right) - \lambda(x_1 - x_2^2)$$
Then we compute the gradient wrt to $x_1,x_2,\lambda$, set to zero and solve.
"""
# ╔═╡ 9a92426e-bc9a-4ead-a93c-ec1c3f23ca57
begin
f(x1,x2) = -exp.(-(x1.*x2 - 3/2).^2 - (x2-3/2).^2)
c(z) = sqrt(z)
x=0:0.01:3.5
end
# ╔═╡ d7acef71-fe0d-4062-9342-9061afde1293
let
p1 = surface(x,x,(x,y)->f(x,y),xlab = L"x_1", ylab = L"x_2")
scatter3d!(p1,[1.358],[1.165],[f(1.358,1.165)],markercolor=:red,leg=false)
p2 = contour(x,x,(x,y)->f(x,y),lw=1.5,levels=[collect(0:-0.1:-0.85)...,-0.887,-0.95,-1],xlab = L"x_1", ylab = L"x_2")
plot!(p2,c,0.01,3.5,label="",lw=2,color=:black)
scatter!(p2,[1.358],[1.165],markersize=5,markercolor=:red,label="Constr. Optimum")
plot(p1,p2,size=(900,300))
end
# ╔═╡ de879e63-fdc9-486e-9aaa-1a55c32d93d5
md"""
* If we had multiple constraints ($l$), we'd just add them up to get
$$
\mathcal{L}(\mathbf{x},\mathbf{\lambda}) = f(\mathbf{x}) - \sum_{i=1}^l \lambda_i h_i(\mathbf{x})$$
"""
# ╔═╡ 54a17a40-fd67-4040-99b9-d0d460385155
md"""
## Inequality Constraints
Suppose now we had
$$\begin{align}
\min_\mathbf{x} & f(\mathbf{x}) \\
\text{subject to } & g(\mathbf{x}) \leq 0
\end{align}$$
which, if the solution lies *on* the constraint *boundary*, means that
$$\nabla f - \mu \nabla g = 0$$
for some scalar $\mu$ - as before.
* In this case, we say the **constraint is active**.
* In the opposite case, i.e. the solution lies **inside** the contrained region, we way the contraint is **inactive**.
* In that case, we are back to an *unconstrained* problem, look for $\nabla f = 0$, and set $\mu=0$.
"""
# ╔═╡ 9f1f0926-b1b3-4ef3-aebe-277531abfd33
let
# the blue area shows the FEASIBLE SET
contour(x,x,(x,y)->f(x,y),lw=1.5,levels=[collect(0:-0.1:-0.85)...,-0.887,-0.95,-1], title = "Constraint is Active or Binding")
plot!(c,0.01,3.5,label="",lw=2,color=:black,fill=(0,0.5,:blue))
scatter!([1.358],[1.165],markersize=5,markercolor=:red,label="Constr. Optimum")
end
# ╔═╡ 2f8d4513-0871-48f3-acfb-f7a0254f7b8e
let
c2(x1) = 1+sqrt(x1)
contour(x,x,(x,y)->f(x,y),lw=1.5,levels=[collect(0:-0.1:-0.85)...,-0.887,-0.95,-1], title = "Constraint is Inactive or Slack")
plot!(c2,0.01,3.5,label="",lw=2,color=:black,fill=(0,0.5,:blue))
scatter!([1],[1.5],markersize=5,markercolor=:red,label="Unconstr. Optimum")
end
# ╔═╡ 8901f1e5-f243-4fe7-afda-6cc951b234a5
md"""
## Infinity Step
* We could do an **infinite step** to avoid *infeasible points*:
$$\begin{align}
f_{\infty\text{-step}} &= \begin{cases}
f(\mathbf{x}) & \text{if } g(\mathbf{x}) \leq 0 \\
\infty & \text{else. }
\end{cases}\\
&= f(\mathbf{x}) + \infty (g(\mathbf{x} > 0)
\end{align}$$
* Unfortunately, this is discontinous and non-differentiable, i.e. hard to handle for algorithms.
* Instead, we use a *linear penalty* $\mu g(\mathbf{x})$ on the objective if the constraint is violated.
"""
# ╔═╡ fe6b7d74-1c6f-4253-9254-57d536bb24c5
md"""
* The penalty provides a lower bound to $\infty$:
$$\mathcal{L}(\mathbf{x},\mu) = f(\mathbf{x}) + \mu g(\mathbf{x})$$
* We can get back the infinite step by maximizing the penalty:
$$f_{\infty\text{-step}} = \max_{\mu\geq 0} \mathcal{L}(\mathbf{x},\mu)$$
* Every infeasible $\mathbf{x}$ returns $\infty$, all others return $f(\mathbf{x})$
"""
# ╔═╡ 19d12886-3980-4116-8772-2c1d17b14270
md"""
## Kuhn-Karush-Tucker (KKT)
* Our problem thus becomes
$$\min_\mathbf{x} \max_{\mu\geq 0} \mathcal{L}(\mathbf{x},\mu)$$
* This is called the **primal problem**. Optimizing this requires:
1. Point is feasible.: $g(\mathbf{x}^*) \leq 0$.
2. Penalty goes into the right direction: $\mu \geq 0$. *Dual feasibility*.
3. Feasible point on the boundary: $\mu g(\mathbf{x}^*) = 0$ has $g(\mathbf{x}) = 0$, otherwise $g(\mathbf{x}) < 0$ and $\mu =0$.
4. Active constraint: $\nabla f(\mathbf{x}^*) - \mu \nabla g(\mathbf{x}^*) = 0$. With an active constraint, we want parallel contours of objective and constraint. When inactive, our optimum just has $\nabla f(\mathbf{x}^*) = 0$, which means $\mu = 0$.
The preceding four conditions are called the **Kuhn-Karush-Tucker (KKT)** conditions. In the above order, and in general terms, they are:
1. Feasibility
2. Dual Feasibility
3. Complementary Slackness
4. Stationarity.
The KKT conditions are the FONCs for problems with smooth constraints.
"""
# ╔═╡ 583a5640-adde-4f4f-a344-4e77e04e431d
md"""
## Duality
We can combine equality and inequality constraints:
$$\mathcal{L}(\mathbf{x},\mathbf{\lambda},\mathbf{\mu}) = f(\mathbf{x}) + \sum_{i} \lambda_i h_i(\mathbf{x}) + \sum_j \mu_j g_j(\mathbf{x})$$
where, notice, we reverted the sign of $\lambda$ since this is unrestricted.
* The Primal problem is identical to the original problem and just as difficult to solve:
$$\min_\mathbf{x} \max_{\mathbf{\mu}\geq 0,\mathbf{\lambda}} \mathcal{L}(\mathbf{x},\mathbf{\mu},\mathbf{\lambda})$$
* The Dual problem reverses min and max:
$$\max_{\mathbf{\mu}\geq 0,\mathbf{\lambda}} \min_\mathbf{x} \mathcal{L}(\mathbf{x},\mathbf{\mu},\mathbf{\lambda})$$
### Dual Values
* The *max-min-inequality* states that for any function $f(a,b)$
$$\max_\mathbf{a} \min_\mathbf{b} f(\mathbf{a},\mathbf{b}) \leq \min_\mathbf{b} \max_\mathbf{a} f(\mathbf{a},\mathbf{b})$$
* Hence, the solution to the dual is a lower bound to the solution of the primal problem.
* The solution to the *dual function*, $\min_\mathbf{x} \mathcal{L}(\mathbf{x},\mathbf{\mu},\mathbf{\lambda})$ is the min of a collection of linear functions, and thus always concave.
* It is easy to optimize this.
* In general, solving the dual is easy whenever minimizing $\mathcal{L}$ wrt $x$ is easy.
"""
# ╔═╡ 9be3cbee-7949-42da-9496-6afc94785f98
md"""
## Penalty Methods
* We can convert the constrained problem back to unconstrained by adding penalty terms for constraint violoations.
* A simple method could just count the number of violations:
$$p_\text{count}(\mathbf{x}) = \sum_{i} (h_i(\mathbf{x}) \neq 0 ) + \sum_j (g_j(\mathbf{x} > 0)$$
* and add this to the objective in an *unconstrained* problem with penalty $\rho > 0$
$$\min_\mathbf{x} f(\mathbf{x}) + \rho p_\text{count}(\mathbf{x})$$
* One can choose the penalty function: for example, a quadratic penaly will produce a smooth objective function
* Notice that $\rho$ needs to become very large sometimes here.
"""
# ╔═╡ 91da671d-7602-490c-a94c-bc9ff4eca59a
md"""
## Interior Point Method
* Also called *barrier method*.
* These methods make sure that the search point remains always feasible.
* As one approaches the constraint boundary, the barrier function goes to infinity. Properties:
1. continuous: $p_\text{barrier}(\mathbf{x})$
2. non negative: $p_\text{barrier}(\mathbf{x})$
3. goes to infinitey : $p_\text{barrier}(\mathbf{x})$ as one approaches the constraint boundary
### Barriers
* Inverse Barrier
$$p_\text{barrier}(\mathbf{x}) = -\sum_i \frac{1}{g_i(\mathbf{x})}$$
* Log Barrier
$$p_\text{barrier}(\mathbf{x}) = -\sum_i \begin{cases}\log(-g_i(\mathbf{x})) & \text{if } g_i(\mathbf{x}) \geq -1 \\
0& \text{else.}
\end{cases}$$
* The approach is as before, one transforms the problem to an unconstrained one and increases $\rho$ until convergence:
$$\min_\mathbf{x} f(\mathbf{x}) + \frac{1}{\rho} p_\text{barrier}(\mathbf{x})$$
### Examples
$$\min_{x \in \mathbb{R}^2} \sqrt{x_2} \text{ subject to }\begin{array}{c} \\
x_2 \geq 0 \\
x_2 \geq (a_1 x_1 + b_1)^3 \\
x_2 \geq (a_2 x_1 + b_2)^3
\end{array}$$
## Constrained Optimisation with [`NLopt.jl`](https://github.com/JuliaOpt/NLopt.jl)
* We need to specify one function for each objective and constraint.
* Both of those functions need to compute the function value (i.e. objective or constraint) *and* it's respective gradient.
* `NLopt` expects contraints **always** to be formulated in the format
$$g(x) \leq 0$$
where $g$ is your constraint function
* The constraint function is formulated for each constraint at $x$. it returns a number (the value of the constraint at $x$), and it fills out the gradient vector, which is the partial derivative of the current constraint wrt $x$.
* There is also the option to have vector valued constraints, see the documentation.
* We set this up as follows:
"""
# ╔═╡ 7e5cedbf-05da-4573-994f-42099fe4b33a
begin
function objective(x::Vector, grad::Vector)
if length(grad) > 0
grad[1] = 0
grad[2] = 0.5/sqrt(x[2])
end
sqrt(x[2])
end
function constraint(x::Vector, grad::Vector, a, b)
if length(grad) > 0
grad[1] = 3a * (a*x[1] + b)^2
grad[2] = -1
end
(a*x[1] + b)^3 - x[2]
end
end
# ╔═╡ b1720238-495e-47f7-ae75-866723f5fca4
x -> x^2
# ╔═╡ 1e31dc48-e849-4045-8dcc-d05889f8609b
let
opt = Opt(:LD_MMA, 2)
lower_bounds!(opt, [-Inf, 0.])
xtol_rel!(opt,1e-4)
min_objective!(opt, objective)
inequality_constraint!(opt, (x,g) -> constraint(x,g,2,0), 1e-8)
inequality_constraint!(opt, (x,g) -> constraint(x,g,-1,1), 1e-8)
(minfunc,minx,ret) = NLopt.optimize(opt, [1.234, 5.678])
end
# ╔═╡ 70893b5f-51e8-4536-b5aa-4e353768bc52
md"""
## NLopt: Rosenbrock
* Let's tackle the rosenbrock example again.
* To make it more interesting, let's add an inequality constraint.
$$\min_{x\in \mathbb{R}^2} (1-x_1)^2 + 100(x_2-x_1^2)^2 \text{ subject to } 0.8 - x_1^2 -x_2^2 \geq 0$$
* in `NLopt` format, the constraint is $x_1^2 + x_2^2 - 0.8 \leq 0$
"""
# ╔═╡ 2a74858c-6d13-4961-86fa-125e4c01e232
function rosenbrockf(x::Vector,grad::Vector)
if length(grad) > 0
grad[1] = -2.0 * (1.0 - x[1]) - 400.0 * (x[2] - x[1]^2) * x[1]
grad[2] = 200.0 * (x[2] - x[1]^2)
end
return (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2
end
# ╔═╡ ec0c21da-c517-4d2e-9832-18eca242621b
function r_constraint(x::Vector, grad::Vector; radius = 0.8)
if length(grad) > 0
grad[1] = 2*x[1]
grad[2] = 2*x[2]
end
return x[1]^2 + x[2]^2 - radius
end
# ╔═╡ ade668a5-87d6-40aa-b670-d0cc41148757
md"""
radius = $(@bind rad Slider(0.1:0.1:2.0,show_value = true,default = 0.2))
"""
# ╔═╡ 52a1ca3b-8a10-4c0a-b4d4-82f995b4b3df
let
grad = zeros(2)
xrange = collect(-2.5:0.01:2)
cc = contour(xrange,xrange, (x,y)->sqrt(rosenbrockf([x, y],grad)), fill=false, color=:viridis, ylab = L"x_2",xlab = L"x_1", leg = :topleft,levels = exp.(range(log(0.1), stop = log(70.0), length = 50)), cbar = false)
contour!(cc,xrange,xrange,(x,y)->r_constraint([x, y],grad,radius = rad), levels = [0],lw = 3)
scatter!(cc,[1],[1],color = :red, lab = "unconstrained optimizer")
# let's compute the constrained optimizer now!
opt = Opt(:LD_MMA, 2)
lower_bounds!(opt, [-5, -5.0])
min_objective!(opt,(x,g) -> rosenbrockf(x,g))
inequality_constraint!(opt, (x,g) -> r_constraint(x,g,radius = rad))
ftol_rel!(opt,1e-9)
(minfunc,minx,ret) = NLopt.optimize(opt, [-1.0,0.0])
scatter!(cc, [minx[1]], [minx[2]], label = "constrained optimizer", leg = :topleft)
end
# ╔═╡ 9cc92292-4531-4f4c-b9b1-2a7a698529ac
md"""
## JuMP.jl
* Introduce [`JuMP.jl`](https://jump.dev/)
* JuMP is a mathematical programming interface for Julia. It is like AMPL, but for free and with a decent programming language.
* The main highlights are:
* It uses automatic differentiation to compute derivatives from your expression.
* It supplies this information, as well as the sparsity structure of the Hessian to your preferred solver.
* It decouples your problem completely from the type of solver you are using. This is great, since you don't have to worry about different solvers having different interfaces.
* In order to achieve this, `JuMP` uses [`MathProgBase.jl`](https://github.com/JuliaOpt/MathProgBase.jl), which converts your problem formulation into a standard representation of an optimization problem.
* Let's look at the readme
* The technical citation is Lubin et al
"""
# ╔═╡ fc75bc7b-0fe0-4f90-a425-6ed6e63acb4a
let
model = Model(GLPK.Optimizer)
@variable(model, 0 <= x <= 2)
@variable(model, 0 <= y <= 30)
# next, we set an objective function
@objective(model, Max, 5x + 3 * y)
# maybe add a constraint called "con"?
@constraint(model, con, 1x + 5y <= 3)
#At this stage JuMP has a mathematical representation of our model internalized
#The MathProgBase machinery knows now exactly how to translate that to different solver interfaces
#For us the only thing left: hit the button!
JuMP.optimize!(model)
(termination_status(model),value(x),value(y))
end
# ╔═╡ 2d7bdf41-85fc-4d0c-a387-60f9c9b36fcb
# JuMP: nonlinear Rosenbrock Example
# Instead of hand-coding first and second derivatives, you only have to give `JuMP` expressions for objective and constraints.
# Here is an example.
let
m = Model(Ipopt.Optimizer)
@variable(m, x)
@variable(m, y)
@NLobjective(m, Min, (1-x)^2 + 100(y-x^2)^2)
@constraint(m, x + y <= 1.5)
JuMP.optimize!(m)
(termination_status(m),value(x),value(y))
end
# ╔═╡ 3ba47370-c4df-4d22-8bda-0b7264927ae6
md"""
## JuMP: Maximium Likelihood
* Let's redo the maximum likelihood example in JuMP.
* Let $\mu,\sigma^2$ be the unknown mean and variance of a random sample generated from the normal distribution.
* Find the maximum likelihood estimator for those parameters!
* density:
$$f(x_i|\mu,\sigma^2) = \frac{1}{\sigma \sqrt{2\pi}} \exp\left(-\frac{(x_i - \mu)^2}{2\sigma^2}\right)$$
* Likelihood Function
$$\begin{aligned}
L(\mu,\sigma^2) = \Pi_{i=1}^N f(x_i|\mu,\sigma^2) =& \frac{1}{(\sigma \sqrt{2\pi})^n} \exp\left(-\frac{1}{2\sigma^2} \sum_{i=1}^N (x_i-\mu)^2 \right) \\
=& \left(\sigma^2 2\pi\right)^{-\frac{n}{2}} \exp\left(-\frac{1}{2\sigma^2} \sum_{i=1}^N (x_i-\mu)^2 \right)
\end{aligned}$$
* Constraints: $\mu\in \mathbb{R},\sigma>0$
* log-likelihood:
$$\log L = l = -\frac{n}{2} \log \left( 2\pi \sigma^2 \right) - \frac{1}{2\sigma^2} \sum_{i=1}^N (x_i-\mu)^2$$
"""
# ╔═╡ 84b81d87-d015-4cc2-8649-c94ba7770bd9
let
distrib = Normal(4.5,3.5)
n = 10000
data = rand(distrib,n);
m = Model(Ipopt.Optimizer)
set_optimizer_attribute(m, MOI.Silent(), true)
@variable(m, μ, start = 0.0)
@variable(m, σ >= 0.0, start = 1.0) # notice bound constraint
@NLobjective(m, Max,
-(n/2) * log(2π * σ^2) - sum( (data[i] - μ)^2 for i = 1:n)/(2 * σ^2)
)
JuMP.optimize!(m)
md"""
parameter | data (truth) | estimate
-------- | ------------ | --------
μ | $(value(μ)) | $(mean(data))
σ | $(value(σ)) | $(std(data))
"""
end
# ╔═╡ 061de138-b6fe-4437-aa26-0afaec2e66f0
md"""
# Linear Constrained Problems (LPs)
* Very similar to before, just that both objective and constraints are *linear*.
$$\begin{align}
\min_\mathbf{x} & \mathbf{c}^T \mathbf{x}\\
\text{subject to } & \mathbf{w}_{LE}^{(i)T} \mathbf{x} \leq b_i \text{ for }i\in{1,2,3,\dots}\\
& \mathbf{w}_{GE}^{(j)T} \mathbf{x} \geq b_j \text{ for }j\in{1,2,3,\dots}\\
& \mathbf{w}_{EQ}^{(k)T} \mathbf{x} = b_k \text{ for }k\in{1,2,3,\dots}\\
\end{align}$$
* Our initial JuMP example was of that sort.
## A Cannery Problem
* A can factory (a cannery) has plants in Seattle and San Diego
* They need to decide how to serve markets New-York, Chicago, Topeka
* Firm wants to
1. minimize shipping costs
2. shipments cannot exceed capacity
3. shipments must satisfy demand
* Formalize that!
* Plant capacity $cap_i$, demands $d_j$ and transport costs from plant $i$ to market $j$, $dist_{i,j} c$ are all given.
* Let $\mathbf{x}$ be a matrix with element $x_{i,j}$ for number of cans shipped from $i$ to $j$.
## From Maths ...
$$\begin{align}
\min_\mathbf{x} & \sum_{i=1}^2 \sum_{j=1}^3 dist_{i,j}c \times x_{i,j}\\
\text{subject to } & \sum_{j=1}^3 x(i,j) \leq cap_i , \forall i \\
& \sum_{i=1}^2 x(i,j) \geq d_j , \forall j
\end{align}$$
## ... to `JuMP`
"""
# ╔═╡ febc5df1-4ddd-48a7-94f8-251d316c72d0
# ... to JuMP
# https://github.com/JuliaOpt/JuMP.jl/blob/release-0.19/examples/cannery.jl
# Copyright 2017, Iain Dunning, Joey Huchette, Miles Lubin, and contributors
# This Source Code Form is subject to the terms of the Mozilla Public
# License, v. 2.0. If a copy of the MPL was not distributed with this
# file, You can obtain one at http://mozilla.org/MPL/2.0/.
#############################################################################
# JuMP
# An algebraic modeling language for Julia
# See http://github.com/JuliaOpt/JuMP.jl
#############################################################################
"""
example_cannery(; verbose = true)
JuMP implementation of the cannery problem from Dantzig, Linear Programming and
Extensions, Princeton University Press, Princeton, NJ, 1963.
Author: Louis Luangkesorn
Date: January 30, 2015
"""
function example_cannery(; verbose = true)
plants = ["Seattle", "San-Diego"]
markets = ["New-York", "Chicago", "Topeka"]
# Capacity and demand in cases.
capacity = [350, 600]
demand = [300, 300, 300]
# Distance in thousand miles.
distance = [2.5 1.7 1.8; 2.5 1.8 1.4]
# Cost per case per thousand miles.
freight = 90
num_plants = length(plants)
num_markets = length(markets)
cannery = Model(GLPK.Optimizer)
@variable(cannery, ship[1:num_plants, 1:num_markets] >= 0)
# Ship no more than plant capacity
@constraint(cannery, capacity_con[i in 1:num_plants],
sum(ship[i,j] for j in 1:num_markets) <= capacity[i]
)
# Ship at least market demand
@constraint(cannery, demand_con[j in 1:num_markets],
sum(ship[i,j] for i in 1:num_plants) >= demand[j]
)
# Minimize transporatation cost
@objective(cannery, Min, sum(distance[i, j] * freight * ship[i, j]
for i in 1:num_plants, j in 1:num_markets)
)
JuMP.optimize!(cannery)
if verbose
println("RESULTS:")
for i in 1:num_plants
for j in 1:num_markets
println(" $(plants[i]) $(markets[j]) = $(JuMP.value(ship[i, j]))")
end
end
end
@test JuMP.termination_status(cannery) == MOI.OPTIMAL
@test JuMP.primal_status(cannery) == MOI.FEASIBLE_POINT
@test JuMP.objective_value(cannery) == 151200.0
end
# ╔═╡ 1334fc66-23bb-4a9a-9534-6f4d7a88d235
example_cannery()
# ╔═╡ 520c3423-01d4-4794-893c-f7e5546d26f6
md"""
# Discrete Optimization / Integer Programming
* Here the choice variable is contrained to come from a discrete set $\mathcal{X}$.
* If this set is $\mathcal{X} = \mathbb{N}$, we have an **integer program**
* If only *some* $x$ have to be discrete, this is a **mixed integer program**
## Example
$$\begin{align}
\min_\mathbf{x} & x_1 + x_2\\
\text{subject to } & ||\mathbf{x}|| \leq 2\\
\text{and } & \mathbf{x} \in \mathbb{N}
\end{align}$$
* continuous optimum is $(-\sqrt{2},-\sqrt{2})$ and objective is $y=-2\sqrt{2}=-2.828$
* Integer constrained problem is only delivering $y=-2$, and $\mathbf{x}^*\in \{(-2,0),(-1,-1),(0,-2)\}$
"""
# ╔═╡ ee374307-043c-473d-b561-192718ea612b
let
x = -3:0.01:3
dx = repeat(range(-3,stop = 3, length = 7),1,7)
contourf(x,x,(x,y)->x+y,color=:heat)
scatter!(dx,dx',legend=false,markercolor=:white)
plot!(x->sqrt(4-x^2),-2,2,c=:white)
plot!(x->-sqrt(4-x^2),-2,2,c=:white)
scatter!([-2,-1,0],[0,-1,-2],c=:red)
scatter!([-sqrt(2)],[-sqrt(2)],c=:red,markershape=:cross,markersize=11)
end
# ╔═╡ 636f7fe1-a518-4f09-bdbd-bcb79703dda5
md"""
## Rounding
* One solution is to just *round the continuous solution to the nearest integer*
* We compute the **relaxed** problem, i.e. the one where $x$ is continuous.
* Then we round up or down.
* Can go terribly wrong.
## Cutting Planes
* This is an exact method
* We solve the relaxed problem first.
* Then we add linear constraints that result in the solution becoming integral.
## Branch and Bound
* This enumerates all possible soultions.
* Branch and bound does this, without having to compute all of them.
"""
# ╔═╡ 2fcf8444-202a-453a-af28-0eade00e3981
md"""
## Example: The Knapsack Problem
* We are packing our knapsack for a trip but only have space for the most valuable items.
* We have $x_i=0$ if item $i$ is not in the sack, 1 else.
$$\begin{align}
\max_x & \sum_{i=1}^n v_i x_i \\
\text{s.t. } & \sum_{i=1}^n w_i x_i \leq w_{max} \\
w_i \in \mathbb{N}_+, & v_i \in \mathbb{R}
\end{align}$$
* If there are $n$ items, we have $2^n$ possible design vectors.
* But there is a useful recursive relationship.
* If we solved $n-1$ knapsack problems so far and deliberate about item $n$
* If it's not worth including item $n$, then the solution **is** the knapsack problem for $n-1$ items and capacity $w_{\max}$
* If it IS worth including it: solution will have value of knapsack with $n-1$ items and reduced capacity, plus the value of the new item
* This is dynamic progamming.
"""
# ╔═╡ 3f231741-12f6-4426-8d46-b6cfe623a248
(1:10).^2
# ╔═╡ 8eef8fc6-2272-46ff-a5ab-b2a4e63d9b7a
let
n , w , v , W = 5, [ 2, 8, 4, 2, 5 ] , [ 5, 3, 2, 7, 4 ] , 10
m = zeros(n,W+1) # m[i,j] is value of including i if capacity is j
for i in 1:n
for (jx,jw) in enumerate(0:W) # if current capacity is
if w[i] > jw
m[i,jx] = i == 1 ? 0.0 : m[i-1,jx] # if weight exceeds current capacity, value of problem is identical to value of problem without the i-th item at that capacity.
else
# weight is feasible to add: should we add it?
m[i,jx] = i == 1 ? v[i] : max(m[i-1,jx], m[i-1, jx - jw] + v[i])
end
end
end
m
end
# ╔═╡ b1f21cfe-ebc2-4321-9636-49739568a7de
md"""
...or, a bit more elegant:
"""
# ╔═╡ b175638b-db07-41c5-a409-f7814d14dc42
let
# Maximization problem
m = Model(GLPK.Optimizer)
set_optimizer_attribute(m, MOI.Silent(), true)
@variable(m, x[1:5], Bin)
profit = [ 5, 3, 2, 7, 4 ]
weight = [ 2, 8, 4, 2, 5 ]
capacity = 10
# Objective: maximize profit
@objective(m, Max, dot(profit, x))
# Constraint: can carry all
@constraint(m, dot(weight, x) <= capacity)
# Solve problem using MIP solver
JuMP.optimize!(m)
OrderedDict( i => Dict(:included => convert(Bool,JuMP.value(x[i])), :profit_over_weight => profit[i]/weight[i]) for i in 1:5 )
end
# ╔═╡ 9644fbd5-298c-4041-822b-1b44917b2c41
begin
danger(head,text) = Markdown.MD(Markdown.Admonition("danger", head, [text]));
danger(text) = Markdown.MD(Markdown.Admonition("danger", "Attention", [text]));
info(text) = Markdown.MD(Markdown.Admonition("info", "Info", [text]));
tip(text) = Markdown.MD(Markdown.Admonition("tip", "Tip", [text]));
midbreak = html"