# JuMP + Ipopt を用いた最適化
- JuMP(Julia for Mathematical Optimization): 数理最適化モデルを記述するためのモデリング言語．簡単だし，ソルバごとにAPIを覚えなくて良い．
- Ipopt(Interior Point OPTimizer): 大規模最適化問題を解くためのオープンソースライブラリ． EPL(Ecliple Public License)で利用できる．主双対内点法を使っているのが特徴．
## JupyterLab と Julia のインストール
- Mac
1. homebrew をインストール
https://brew.sh/ を参考にすればOK
2. JupyterLab と Julia をインストール
homebrew が入っている状態でターミナルから以下を実行：

       brew install --cask jupyterlab
       brew install --cask julia

3. Julia から JupyterLab を使えるようにする
ターミナルから Julia を起動し， `]`をタイプしてパッケージモードにした後,

        activate .
        add IJulia
        build IJulia

を実行． [ここ](https://qiita.com/shyu_manabe/items/3978c1ef5d96d4e9dcef)とか[
ここ](https://julia.quantecon.org/getting_started_julia/getting_started.html#installing-packages)が参考になる． `julia 1.10`からは上のコマンドを使わないとうまく動かないようだ．
- Windows
Julia と Jupyterlab をそれぞれインストール．
たとえば[ここ](https://zenn.dev/ohno/books/356315a0e6437c/viewer/c7687c)とかを参考にする
## JuMPとIpOptのインストール
JupyterLab のセルに以下を記述しておく．

import Pkg
Pkg.add("JuMP")
Pkg.add("Ipopt")

# 最適化をやってみる
## 最初のサンプル

using JuMP, Ipopt # モジュールを使用
model = Model(Ipopt.Optimizer) # モデルを定義
@variable(model, x >= 0) # 変数xを追加
@variable(model, 0 <= y <= 3) # 変数 yを追加
@objective(model, Min, 12x + 20y) # 目的関数を設定
@constraint(model, c1, 6x + 8y >= 100) # 1つめの制約条件を追加
@constraint(model, c2, 7x + 12y >= 120) # 2つめの制約条件を追加
print(model) # 構築されたモデルを数式表示
optimize!(model) # 最適化

これだけで最適化される． 下記を実行すれば計算結果を出力させられる．

display( termination_status(model) )
display( primal_status(model) )
display( dual_status(model) )
display( objective_value(model) )
display( value(x) )
display( value(y) )
display( shadow_price(c1) )
display( shadow_price(c2) )

## 利用者均衡配分
$K$本の平行リンクで構成されるネットワークを考え,起点から終点への総交通需要を$Q$で表す．リンク交通量を$\boldsymbol{x} = (x_{k}: k =1, \cdots, K)$で表し，各リンクの所要時間が線形式：

$$
\boldsymbol{t}(\boldsymbol{x}) = \boldsymbol{D} \boldsymbol{x} + \boldsymbol{e}
$$

で表されるとする．ここで，$\boldsymbol{D}$は正定値対称な$K\times{}K$行列とし，$\boldsymbol{e}$は$K$次元列ベクトルである．このとき，利用者均衡配分問題は，以下の最適化問題に帰着する：

\begin{align}
\min_{\boldsymbol{x}} &f(\boldsymbol{x}) = \frac{1}{2}\boldsymbol{x}^{\top} \boldsymbol{D} \boldsymbol{x} + \boldsymbol{e}^{\top}\boldsymbol{x}\\
\text{Sub to.} \quad & \boldsymbol{1}^{\top} \boldsymbol{x}= Q\\
& \boldsymbol{x} \geq \boldsymbol{0}
\end{align}

using Random # Dとeをランダムに生成するために Random と LinearAlgebraを使う
using LinearAlgebra
K = Int(5e4) # 経路数
Q = 1 # 総交通量
D = Diagonal([rand(1:0.1:10) for k in 1:K]) # リンク所要時間係数
e = Vector([rand(1:0.1:10) for k in 1:K]) # リンク所要時間定数
f(x) = 0.5x'*D*x + x'*e # 目的関数
∇f(x) = D*x+e # 目的関数の勾配
∇²f(x) = D # 目的関数の Hessian
model = Model(Ipopt.Optimizer) # モデルを定義
set_attribute(model, "print_level", 3) # 出力レベルを設定
@variable(model, x[1:K] .>= 0) # K次元の非負未知変数を定義
@constraint(model, c1, sum(x) == Q) # 総交通量保存則
@operator(model, op_f, K, f, ∇f, ∇²f) # 目的関数fに勾配∇fと Hessian ∇2f を明示的に指定した演算 op_f を定義
@objective(model, Min, op_f(x)) # 目的関数に op_f を指定
optimize!(model) # 最適化
_x = value.(x) # 計算結果
_rho = -shadow_price(c1)
_obj = objective_value(model)
sum(_x.*(∇f(_x).-_rho))