using JuMP
using Ipopt
using LinearAlgebra
"""
Problem: Solve the Rosenbrock function
"""
# For more information see here:
# https://jump.readthedocs.org/en/latest/installation.html#getting-solvers
# Information about Ipopt options here:
# http://www.coin-or.org/Ipopt/documentation/node39.html
m = Model(with_optimizer(Ipopt.Optimizer))
# setsolver(m, IpoptSolver(tol = 1e-6, max_iter = 200, output_file = "results.txt"))
@variable(m, x1)
@variable(m, x2)
@NLobjective(m, Min, 100(x2-x1^2)^2 + (1-x1)^2)
set_start_value(x1, -1.2)
set_start_value(x2, 1)
optimize!(m)
println("got ", objective_value(m), " at [", value(x1), ", ", value(x2), "]")
print(m)
m = Model(with_optimizer(Ipopt.Optimizer))
@variable(m, x[i=1:2])
@NLobjective(m, Min, 100(x[2]-x[1]^2)^2 + (1-x[1])^2)
set_start_value(x[1], -1.2)
set_start_value(x[2], 1)
println(m)
optimize!(m)
println("got ", objective_value(m), " at ", [value(x[1]), value(x[2])])
include("mle.jl")
n = 100
m = Model(with_optimizer(Ipopt.Optimizer))
@variable(m, x[i=1:n], start=(i/2))
@NLexpression(m, myexpr[i=1:n], sin(x[i]))
@NLconstraint(m, myconstr[i=1:n], myexpr[i] <= 0.5)
@NLobjective(m, Min, sum(x[i]^2-1 for i = 1:n))
m
optimize!(m)
objective_value(m)
sol = [value(x[i]) for i =1:100]
sin.(sol)
n = 100
m = Model(with_optimizer(Ipopt.Optimizer))
@variable(m, x[i=1:n], start=-1.0)
@NLexpression(m, myexpr[i=1:n], sin(x[i]))
@NLconstraint(m, myconstr[i=1:n], myexpr[i] <= 0.5)
@NLobjective(m, Min, sum(x[i]^2-1 for i = 1:n))
optimize!(m)
sol = [value(x[i]) for i =1:100]
objective_value(m)
sin.(sol)
Consider the program \begin{align*} \min_x\ & c^Tx \\ \mbox{s.t. } & Ax = b \\ & x \geq 0 \end{align*}
The Lagrangian is $$ L(x,\lambda,\mu) = c^Tx + \lambda^T(Ax - b) + \mu^T(-x) $$ and the duality bound is given by $$ \min_{x \in S} \max_{\mu \geq 0, \lambda} L(x,\lambda, \mu) $$ The KKT system becomes \begin{align*} c + A^T\lambda - \mu &= 0 \\ Ax &= b \\ x &\geq 0 \\ \mu_i x_i &= 0,\ i = 1,\ldots,n \\ \mu &\geq 0 \end{align*} Considering $\mu$ as a surplus vector, the first and fifth conditions can be rewritten as $$ c + A^T\lambda \geq 0 $$
At $(x^*, \lambda^*, \mu^*)$, $$ L(x^*, \lambda^*, \mu^*) = c^Tx^* + (\lambda^*)^T(Ax^*-b)-(\mu^*)^Tx^* = c^Tx^* $$
Moreover, the first KKT condition implies that at the primal-dual solution $(x^*, \lambda^*, \mu^*)$, the Lagrangian becomes $$ L(x^*, \lambda^*, \mu^*) = -b^T\lambda^* $$ as \begin{align*} L(x^*, \lambda^*, \mu^*) = c^Tx^* + (\lambda^*)^T(Ax^*-b)-(\mu^*)^Tx^* & = (c^T+(\lambda^*)^TA-(\mu^*)^T)x^*-(\lambda^*)^Tb \\ & = (x^*)^T(c+A^T\lambda^*-\mu^*)-b^T\lambda^* \end{align*}
Along with the duality bound, this means that we have to solve the problem \begin{align*} \max_{\lambda} \ & -b^T\lambda \\ \mbox{s.t. } & A^T\lambda \geq -c \end{align*} Setting $y = -\lambda$, the program can be rewritten as \begin{align*} \max_{y}\ & b^Ty \\ \mbox{s.t. } & A^Ty \leq c \end{align*} We retrieve the well-known dual problem is linear programming.
If we apply the logarithmic barrier to the primal linear program, we can rewrite it as \begin{align*} \min\ & c^Tx - \mu \sum_{i = 1}^n \ln x_i \\ \mbox{s.t. } & Ax = b, \end{align*} and the barrier Lagrangian is $$ L(x,\lambda) = c^Tx + \lambda^T(Ax - b) - \mu \sum_{i = 1}^n \ln x_i $$ and the KKT conditions are now \begin{align*} c + A^T\lambda - \mu X^{-1}e &= 0 \\ Ax &= b \end{align*} Set $s = \mu X^{-1}e$. We can rewrite the conditions as \begin{align*} c + A^T\lambda - s &= 0 \\ Ax &= b \\ s &= \mu X^{-1}e \end{align*} or \begin{align*} c + A^T\lambda - s &= 0 \\ Ax &= b \\ s_i x_i &= \mu \end{align*} The logarithmic barrier is equivalent to perturb the complementarity condition of the KKT system.
As before, set $y = -\lambda$. Solving the KKT system is equivalent to look for a solution of the linear system \begin{align*} A^Ty+s-c &= 0\\ Ax-b &= 0 \\ SXe-\mu e &= 0 \end{align*} The corresponding Newton equation is then $$ \begin{pmatrix} 0 & A^T & I \\ A & 0 & 0 \\ S & 0 & X \end{pmatrix} \begin{pmatrix} \Delta x \ \Delta y \ \Delta s
- \begin{pmatrix} A^Ty+s-c \\ Ax-b \\ SXe-\mu e \end{pmatrix} $$
Adapted from https://www.sonoma.edu/users/w/wilsonst/Courses/Math_131/lp/Farm.html
A farmer has 10 acres to plant in wheat and rye. He has to plant at least 7 acres. However, he has only CAD1200 to spend and each acre of wheat costs CAD200 to plant and each acre of rye costs CAD100 to plant. Moreover, the farmer has to get the planting done in 12 hours and it takes an hour to plant an acre of wheat and 2 hours to plant an acre of rye. If the profit is CAD500 per acre of wheat and CAD300 per acre of rye how many acres of each should be planted to maximize profits?
c = [-500.0 ; -300.0 ; 0.0 ; 0.0 ; 0.0; 0.0]
A = [1.0 1.0 1.0 0.0 0.0 0.0; 1.0 1.0 0.0 -1.0 0.0 0.0; 2.0 1.0 0.0 0.0 1.0 0.0; 1.0 2.0 0.0 0.0 0.0 1.0]
b = [10.0; 7.0; 12.0; 12]
x = [3.5; 4; 2.5; 0.5; 1; 0.5]
m, n = size(A)
μ = 10.0
y = zeros(4)
s = zeros(6)
for i = 1:6
s[i] = μ/x[i]
end
w = [ x ; y; s ]
θ = 0.5
ρ = 0.9
#using Gurobi
using Clp
model = Model(with_optimizer(Clp.Optimizer))
n = length(x)
m = length(b)
@variable(model, X[1:n] >= 0)
for i = 1:m
@constraint(model, sum(A[i,j]*X[j] for j in 1:n) == b[i])
end
@objective(model, Min, sum(c[i]*X[i] for i in 1:n))
model
optimize!(model)
objective_value(model)
[value(X[i]) for i = 1:n]
cost(c::Vector,x::Vector) = c'*x
cost(c,x)
A*x-b
s, x, s.*x
m,n = size(A)
B = [ zeros(n,n) A' I ; A zeros(m,m) zeros(m,n) ; diagm(0 => x) zeros(n,m) diagm(0 => s)]
This is absolutely inefficient as the matrix is very sparse!
rhs = [ c-s-A'*y ; b-A*x ; [μ - x[i]*s[i] for i=1:n] ]
function d(c::Vector, A::Matrix, b::Vector, x::Vector, y::Vector, s::Vector, μ::Float64)
m,n = size(A)
B = [ zeros(n,n) A' I ; A zeros(m,m) zeros(m,n) ; diagm(0 => x) zeros(n,m) diagm(0 => s)]
rhs = [ c-s-A'*y ; b-A*x ; [μ - x[i]*s[i] for i=1:n] ]
Δ = B\rhs
return Δ
end
Δ = d(c,A,b,x,y,s,μ)
Δx = Δ[1:n]
Δy = Δ[n+1:n+m]
Δs = Δ[n+m+1:n+m+n]
A*(x+Δx)
x+Δx, cost(c,x+Δx)
function maxα(x::Vector, s::Vector, Δx::Vector, Δs::Vector)
n = length(c)
α = Inf
for i=1:n
if (Δx[i] < 0)
α = min(α,-x[i]/Δx[i])
end
if (Δs[i] < 0)
α = min(α,-s[i]/Δs[i])
end
end
return α
end
α = maxα(x,s,Δx,Δs)
x+α*Δx
cost(c,x+α*Δx)
s+α*Δs
α *= θ
x += α*Δx
y += α*Δy
s += α*Δs
x, y, s
A*x
cost(c,x)
μ *= ρ
while (μ > 1e-50)
Δ = d(c,A,b,x,y,s,μ)
Δx = Δ[1:n]
Δy = Δ[n+1:n+m]
Δs = Δ[n+m+1:n+m+n]
α = θ*maxα(x,s,Δx,Δs)
x += α*Δx
y += α*Δy
s += α*Δs
μ *= ρ
end
x
cost(c,x)