{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Juliaによる混合整数計画法（線形計画法）\n",
"\n",
"最近~~仕事さぼって~~Juliaでいろんな実験をするのにはまっています。今日はJuliaと外部APIを使って混合整数計画法をやってみます。\n",
"\n",
"私は線形計画法のプロでもなんでもないので，線形計画法については[こちら（PDF直リンク）](http://web.tuat.ac.jp/~miya/fujie_ORSJ.pdf)が参考になると思います。[Rによる項目反応理論](https://shop.ohmsha.co.jp/shopdetail/000000003874/)の第15章にも取り上げられている内容です。\n",
"\n",
"線形計画法とは特定の制約条件(Constraints)下で目的関数(Objective)を最小or最大とするような変数の組み合わせを求めるもので，\n",
"\n",
"$$\n",
"\\begin{alignat}{2}\n",
" &\\text{minimize} & \\quad c^\\top x & \\label{eqn:P_Obj} \\\\\n",
" &\\text{subject to} & \\quad Ax &\\geq b \\label{eqn:P_Con-Eq} \\\\\n",
" & & \\quad x &\\geq 0. \\label{eqn:P_Con-Non} \n",
"\\end{alignat}\n",
"$$\n",
"\n",
"という感じで定式化されます。単純な線形計画法は中学か高校の数学でも習いますが，今回はこれに整数条件が加わった整数計画法および混合整数計画法を使って，「理想的な自動テスト編成」をやります。"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 項目反応理論とテスト情報量\n",
"\n",
"項目反応理論（IRT)ではテスト項目やテスト自体の測定精度の指標として**テスト情報量 (Test Information, TIF)** が定義されています。TIFはICC（項目特性曲線）を $ P(\\theta) = logistic(\\alpha(\\theta-\\beta)) $ とおき，あるテストに含まれるテスト項目に関する添え字を $j = 1, 2, ... J$とすると，以下のように求められます。\n",
"$$\n",
"TIF(\\theta) = \\sum_j^J \\alpha_j^2(p(\\theta)(1-P(\\theta))\n",
"$$\n",
"\n",
"TIFは受検者の潜在能力値 $\\theta$ の関数として定義され，テスト情報量が大きいところほど，測定精度が高いと解釈されます（ちなみにテスト情報量自体は統計学におけるいわゆる**フィッシャー情報量**の定義から導かれます。対数尤度関数の二階偏微分の負の期待値ですね）。\n",
"\n",
"ここで，あるテスト項目プール（IRTのパラメタ推定済みの項目をたくさん保持しておく場所，項目バンクとも）から，一定の条件でテスト冊子を編成するケースを考えます。たとえば100個の問題の中から，出題する対象に最適な測定精度のテスト冊子を構成するのです。\n",
"\n",
"最適化の制約や方向には様々な方法があるのですが，ここではテスト情報量の**下限値制約**を満たすなかで，最もテスト項目の**情報量を節約**しつつ，**重複項目のないテスト冊子**を作ることを目的にしたいと思います。下限値制約は一定以上の測定精度を保持するため，情報量の節約は項目バンクの暴露（IRTでは項目が受検者に対策されてしまうと，推定が正しくおこなえなくなる恐れがある）を防ぐためです。情報量の節約は無駄な項目の出題を避けることもできるので，効率よく受検者の能力を推定することにも役に立ちます。"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Solver\n",
"\n",
"基本的にはJuliaの`JuMP`pkgを使います。しかし最適化のsolverに関してはいくつか選択肢があり，非商用から商用まで様々カバーしています。\n",
"\n",
"ここでは，変な牛さんみたいなアイコンが特徴の[GLPK](http://www.gnu.org/software/glpk/)を使います。"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"## DO NOT RUN ## \n",
"# using Pkg\n",
"# Pkg.add(\"JuMP\"); Pkg.add(\"GLPK\")"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"tif (generic function with 1 method)"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using JuMP, GLPK\n",
"using Random, Distributions, StatsFuns # 乱数生成と情報量計算に使う\n",
"using Plots\n",
"\n",
"# test infromation function\n",
"function tif(θ, α, β)\n",
" I = zeros(length(θ))\n",
" for m in 1:length(θ)\n",
" for j in 1:length(α)\n",
" p = logistic(α[j]*(θ[m]-β[j]))\n",
" I[m] += α[j]^2 * p * (1-p)\n",
" end\n",
" end\n",
" return I\n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 自動テスト編成 その1\n",
"\n",
"まず手始めに，90個の項目をもつ項目プールから，30個のテスト冊子を取り出すことを考えます。最低限度のテスト情報量を保持したテスト冊子を目指します。目的となるテスト情報量は項目バンク全体のテスト情報量を，取り出すテスト冊子の項目数に合わせて除したものとします。つまり，もとの項目バンクとなるべく同じテスト情報関数の形状になるように，テスト冊子を構成するのです。詳細な制約条件はまた後ほど記述します。"
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 34,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"Random.seed!(1234)\n",
"α = rand(LogNormal(-0.5, 0.5), 90) # 項目識別力\n",
"histogram(α)"
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 35,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"β = rand(Normal(0, 1.5), 90) # 項目困難度\n",
"scatter(α, β)"
]
},
{
"cell_type": "code",
"execution_count": 48,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 48,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"target_θ = [-2:1:2;]\n",
"target_TIF = tif(target_θ, α, β) * 30/90\n",
"IIF = [tif(target_θ[i], α[j], β[j])[1] for i in 1:5, j in 1:90]\n",
"plot([-4:0.1:4;], tif([-4:0.1:4;], α, β), label = \"target TIF\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"最適化方向に対しての情報量の下限値は適当な点数でOKです。このtarget_TIFが変な形（きつい三角形）だったりすると最適化が難しく，失敗することがあります。\n",
"\n",
"なお，最適化のための目的関数ですが，項目バンクのテスト情報量をなるべく節約する方向に設定すれば良いので，**テスト冊子に含まれる識別力の和**とすればよいです。もちろん無制約に最小化すれば，ただ識別力の小さい項目をとってくるだけですので，ここに下限値制約が加わるわけです。"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"続いて，モデルを組み立てていきます。必要になるのは\n",
"1. 目的関数\n",
"2. 決定変数\n",
"3. 制約条件\n",
"の3つです。\n",
"\n",
"`JuMP`の文法については[公式ドキュメント](http://www.juliaopt.org/JuMP.jl/v0.19.2/)やsolverへのインターフェイスを提供する関数の[ドキュメント](http://www.juliaopt.org/MathOptInterface.jl/v0.6.1/apimanual.html#Sets-and-Constraints-1)が参考になります。日本語だと[この記事](https://myenigma.hatenablog.com/entry/2017/10/28/141805)も参考になります。"
]
},
{
"cell_type": "code",
"execution_count": 38,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$$ 0.9358253619513699 x_{1} + 0.3864039678519896 x_{2} + 0.47367238044788285 x_{3} + 0.3861778953382992 x_{4} + 0.9344479579417156 x_{5} + 1.832972103656246 x_{6} + 0.7916836815603108 x_{7} + 0.5294758637815399 x_{8} + 0.7797103676313849 x_{9} + 0.46837228882930027 x_{10} + 0.458291117488107 x_{11} + 0.6007082574478093 x_{12} + 0.6466385585208436 x_{13} + 1.5317202938948689 x_{14} + 0.40096476799044956 x_{15} + 0.6408550768129122 x_{16} + 0.5349468568503699 x_{17} + 0.7296845347024683 x_{18} + 0.628800134841628 x_{19} + 0.28601393983292045 x_{20} + 1.3258902462243427 x_{21} + 0.3016862576970978 x_{22} + 1.0541649316862556 x_{23} + 0.3487622026564447 x_{24} + 0.12176286868804573 x_{25} + 0.5844948721208734 x_{26} + 0.6540887775920585 x_{27} + 0.8910445348777064 x_{28} + 0.5194024376210554 x_{29} + 0.4487212333196793 x_{30} + 0.3198714428156127 x_{31} + 0.9986594771333122 x_{32} + 0.7055421634995673 x_{33} + 0.5955779518373137 x_{34} + 0.6511515482808155 x_{35} + 0.7871285738761749 x_{36} + 0.9496838982854913 x_{37} + 0.46918180201107584 x_{38} + 0.4137889356458933 x_{39} + 0.28063140632430184 x_{40} + 0.5827009100968081 x_{41} + 0.35147749087145375 x_{42} + 0.4537275947780989 x_{43} + 0.5180318132597422 x_{44} + 0.3070564068919007 x_{45} + 0.5727943082397413 x_{46} + 0.6589672012247402 x_{47} + 0.49449461582580767 x_{48} + 0.3660842578939421 x_{49} + 0.46213305833115575 x_{50} + 0.32845272504237755 x_{51} + 0.4626160382997182 x_{52} + 0.43031113238387964 x_{53} + 0.4246601908750506 x_{54} + 0.5150303499858949 x_{55} + 0.7845994732834861 x_{56} + 2.031423286906088 x_{57} + 0.519968612413358 x_{58} + 1.1304905944970691 x_{59} + 0.5915700850197081 x_{60} + 0.35794035521067413 x_{61} + 0.3702457117451199 x_{62} + 1.2726307109702621 x_{63} + 0.4670186552738027 x_{64} + 0.2751743485977304 x_{65} + 0.6478611794323723 x_{66} + 0.7585693747069833 x_{67} + 0.49752705480509524 x_{68} + 0.7286124025636858 x_{69} + 0.8276512803754268 x_{70} + 0.6645095670616384 x_{71} + 1.701935676118738 x_{72} + 0.2990142496524239 x_{73} + 0.6487146800256491 x_{74} + 0.4167743821707079 x_{75} + 0.6605971852632762 x_{76} + 0.2745942263865896 x_{77} + 0.2259367360235037 x_{78} + 0.6301540154430951 x_{79} + 0.9045366846224644 x_{80} + 0.4201763847878517 x_{81} + 0.6906205879234821 x_{82} + 0.3590457199557009 x_{83} + 0.6701890578955982 x_{84} + 0.4051936235936281 x_{85} + 0.669584486975328 x_{86} + 0.3403899642463468 x_{87} + 0.913341832689808 x_{88} + 0.6541981558411987 x_{89} + 1.1315811374546587 x_{90} $$"
],
"text/plain": [
"0.9358253619513699 x[1] + 0.3864039678519896 x[2] + 0.47367238044788285 x[3] + 0.3861778953382992 x[4] + 0.9344479579417156 x[5] + 1.832972103656246 x[6] + 0.7916836815603108 x[7] + 0.5294758637815399 x[8] + 0.7797103676313849 x[9] + 0.46837228882930027 x[10] + 0.458291117488107 x[11] + 0.6007082574478093 x[12] + 0.6466385585208436 x[13] + 1.5317202938948689 x[14] + 0.40096476799044956 x[15] + 0.6408550768129122 x[16] + 0.5349468568503699 x[17] + 0.7296845347024683 x[18] + 0.628800134841628 x[19] + 0.28601393983292045 x[20] + 1.3258902462243427 x[21] + 0.3016862576970978 x[22] + 1.0541649316862556 x[23] + 0.3487622026564447 x[24] + 0.12176286868804573 x[25] + 0.5844948721208734 x[26] + 0.6540887775920585 x[27] + 0.8910445348777064 x[28] + 0.5194024376210554 x[29] + 0.4487212333196793 x[30] + 0.3198714428156127 x[31] + 0.9986594771333122 x[32] + 0.7055421634995673 x[33] + 0.5955779518373137 x[34] + 0.6511515482808155 x[35] + 0.7871285738761749 x[36] + 0.9496838982854913 x[37] + 0.46918180201107584 x[38] + 0.4137889356458933 x[39] + 0.28063140632430184 x[40] + 0.5827009100968081 x[41] + 0.35147749087145375 x[42] + 0.4537275947780989 x[43] + 0.5180318132597422 x[44] + 0.3070564068919007 x[45] + 0.5727943082397413 x[46] + 0.6589672012247402 x[47] + 0.49449461582580767 x[48] + 0.3660842578939421 x[49] + 0.46213305833115575 x[50] + 0.32845272504237755 x[51] + 0.4626160382997182 x[52] + 0.43031113238387964 x[53] + 0.4246601908750506 x[54] + 0.5150303499858949 x[55] + 0.7845994732834861 x[56] + 2.031423286906088 x[57] + 0.519968612413358 x[58] + 1.1304905944970691 x[59] + 0.5915700850197081 x[60] + 0.35794035521067413 x[61] + 0.3702457117451199 x[62] + 1.2726307109702621 x[63] + 0.4670186552738027 x[64] + 0.2751743485977304 x[65] + 0.6478611794323723 x[66] + 0.7585693747069833 x[67] + 0.49752705480509524 x[68] + 0.7286124025636858 x[69] + 0.8276512803754268 x[70] + 0.6645095670616384 x[71] + 1.701935676118738 x[72] + 0.2990142496524239 x[73] + 0.6487146800256491 x[74] + 0.4167743821707079 x[75] + 0.6605971852632762 x[76] + 0.2745942263865896 x[77] + 0.2259367360235037 x[78] + 0.6301540154430951 x[79] + 0.9045366846224644 x[80] + 0.4201763847878517 x[81] + 0.6906205879234821 x[82] + 0.3590457199557009 x[83] + 0.6701890578955982 x[84] + 0.4051936235936281 x[85] + 0.669584486975328 x[86] + 0.3403899642463468 x[87] + 0.913341832689808 x[88] + 0.6541981558411987 x[89] + 1.1315811374546587 x[90]"
]
},
"execution_count": 38,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"mod1 = Model(with_optimizer(GLPK.Optimizer))\n",
"@variable(mod1, 0 <= x[1:90] <= 1, Int) # 決定変数\n",
"@objective(mod1, Min, α' * x) # 目的変数"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"本来決定変数は0か1の二値です（プールの項目がテスト冊子に含まれるのならば1，そうでないならば0）。`@variavle`の三つ目の引数で整数制約を課しています。本当なら`Int`ではなく`Bin`でもいけそうなのですが，私の環境では無理でした。\n",
"\n",
"当然ですが，目的変数に決定変数が入るわけですから，決定変数は目的変数よりも先に宣言しておく必要があります。\n",
"\n",
"出力されたベクトルが，目的関数の数式表現です。"
]
},
{
"cell_type": "code",
"execution_count": 41,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$ x_{30} + x_{60} + x_{90} \\in \\[0.0, 1.0\\] $"
],
"text/plain": [
"x[30] + x[60] + x[90] in [0.0, 1.0]"
]
},
"execution_count": 41,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"@constraint(mod1, fill(1, 90)'x == 30) # 項目数は30個\n",
"@constraint(mod1, [fill(1,50);fill(0,40)]'x == 18) # 前半から18項目\n",
"@constraint(mod1, [fill(0,50);fill(1,40)]'x == 12) # 後半から12項目\n",
"@constraint(mod1, target_TIF[1]+.1 >= IIF[1,:]'x >= target_TIF[1]-0) # テスト情報量の制約条件1~5\n",
"@constraint(mod1, target_TIF[2]+.1 >= IIF[2,:]'x >= target_TIF[2]-0)\n",
"@constraint(mod1, target_TIF[3]+.1 >= IIF[3,:]'x >= target_TIF[3]-0)\n",
"@constraint(mod1, target_TIF[4]+.1 >= IIF[4,:]'x >= target_TIF[4]-0)\n",
"@constraint(mod1, target_TIF[5]+.1 >= IIF[5,:]'x >= target_TIF[5]-0)\n",
"@constraint(mod1, 1 >= [fill(0,29);1;fill(0,29);1;fill(0, 29);1]'x >= 0) # 敵対項目"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"制約条件は以下の通りです。\n",
"\n",
"- テスト項目数は30とする。\n",
"- 項目バンクの前半と後半を異なる領域の問題（たとえば基礎と応用）とみなして，前半から18項目，後半から12項目採用する。\n",
"- 5点のテスト情報量において，目的情報量よりも+0.1の範囲の誤差に抑える。\n",
"- 敵対項目（一緒に含めてはいけない項目）が3項目（30，60，90）あるとする。\n",
"\n",
"for loopを使えばベタうちしなくても良いと思います。\n",
"\n",
"それでは最適化を実行します。"
]
},
{
"cell_type": "code",
"execution_count": 42,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 0.209900 seconds (23.15 k allocations: 1.387 MiB)\n"
]
}
],
"source": [
"@time optimize!(mod1)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"コンパイル済みならば瞬殺です。結果がうまくいっているかは`termination_status()`で確認できます"
]
},
{
"cell_type": "code",
"execution_count": 43,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"OPTIMAL::TerminationStatusCode = 1"
]
},
"execution_count": 43,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"termination_status(mod1) # OPTIMAL::TerminationStatusCode = 1 ならいちおうOK"
]
},
{
"cell_type": "code",
"execution_count": 47,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 47,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"get_item = value.(x)\n",
"plot([-4:0.1:4;], tif([-4:0.1:4;], α[isone.(get_item)], β[isone.(get_item)]), label = \"optimized test\")\n",
"plot!([-4:0.1:4;], tif([-4:0.1:4;], α, β) * 30/90,label = \"target TIF\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"最適化で指定した点の位置(-2, -1, 0, 1, 2)では良い感じです。端っこの方の領域でズレが生じていますが，どのみにここら辺の帯域の受検者は今回のテストの想定する対象外なので，問題ないです。\n",
"\n",
"### 重複を含まずに，複数のテストフォームを取り出す\n",
"\n",
"次はもっと複雑な条件にします。1つの項目プールから，3つのテストフォームを同時に編成することを考えます。もちろん先ほどと同様にテスト情報量に関する制約があります。さらにフォーム間で項目が重複してはいけないものとします。\n",
"\n",
"さっそくモデルを組み立てていきます。"
]
},
{
"cell_type": "code",
"execution_count": 52,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 52,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# 項目パラメタの生成\n",
"Random.seed!(1234)\n",
"α2 = rand(LogNormal(-0.5, 0.6), 60)\n",
"β2 = rand(Normal(0, 1.5), 60)\n",
"# 制約のためのIIF, TIFの計算\n",
"IIF3 = [tif(target_θ[i], α2[j], β2[j])[1] for i in 1:5, j in 1:60]\n",
"target_TIF3 = tif(target_θ, α2, β2) * 10/60\n",
"plot([-4:0.1:4;], tif([-4:0.1:4;], α2, β2), label = \"target TIF\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"同時複数テスト編成になると，組み合わせの数が一気に跳ね上がり，計算時間が途方もないことになります。そのため今回は全体的に小規模なテストを想定しています。"
]
},
{
"cell_type": "code",
"execution_count": 78,
"metadata": {},
"outputs": [],
"source": [
"# create model\n",
"mod2 = Model(with_optimizer(GLPK.Optimizer))\n",
"# declare variables\n",
"@variable(mod2, 0 <= F1[1:60*3] <= 1, Int) # Test From 1 to 3\n",
"@variable(mod2, 0 <= F2[1:60*3] <= 1, Int)\n",
"@variable(mod2, 0 <= F3[1:60*3] <= 1, Int)\n",
"# objective functions\n",
"@objective(mod2, Min, repeat(α2, inner = 1, outer = 3)' * F1)\n",
"@objective(mod2, Min, repeat(α2, inner = 1, outer = 3)' * F2)\n",
"@objective(mod2, Min, repeat(α2, inner = 1, outer = 3)' * F3)\n",
"# impose constraints\n",
"# item size in each test forms\n",
"@constraint(mod2, [fill(1,60); fill(0,120) ]'F1 == 10) # item size\n",
"@constraint(mod2, [fill(0,60); fill(1,60); fill(0,60)]'F2 == 10)\n",
"@constraint(mod2, [fill(0,120); fill(1, 60)]'F3 == 10)\n",
"# item domain\n",
"@constraint(mod2, [fill(1,35);fill(0,25);fill(0,120)]'F1 == 7) # 前半から18項目\n",
"@constraint(mod2, [fill(0,35);fill(1,25);fill(0,120)]'F1 == 3) # 後半から12項目\n",
"@constraint(mod2, [fill(0,60);fill(1,35);fill(0,25);fill(0,60)]'F2 == 7)\n",
"@constraint(mod2, [fill(0,60);fill(0,35);fill(1,25);fill(0,60)]'F2 == 3)\n",
"@constraint(mod2, [fill(0,120);fill(1,35);fill(0,25)]'F3 == 7)\n",
"@constraint(mod2, [fill(0,120);fill(0,35);fill(1,25)]'F3 == 3)\n",
"# TIF constraints\n",
"@constraint(mod2, target_TIF3[1]+.1 >= [IIF3[1,:];fill(0,120)]'F1 >= target_TIF3[1]-.1)\n",
"@constraint(mod2, target_TIF3[2]+.1 >= [IIF3[2,:];fill(0,120)]'F1 >= target_TIF3[2]-.1)\n",
"@constraint(mod2, target_TIF3[3]+.1 >= [IIF3[3,:];fill(0,120)]'F1 >= target_TIF3[3]-.1)\n",
"@constraint(mod2, target_TIF3[4]+.1 >= [IIF3[4,:];fill(0,120)]'F1 >= target_TIF3[4]-.1)\n",
"@constraint(mod2, target_TIF3[5]+.1 >= [IIF3[5,:];fill(0,120)]'F1 >= target_TIF3[5]-.1)\n",
"\n",
"@constraint(mod2, target_TIF3[1]+.1 >= [fill(0,60);IIF3[1,:];fill(0,60)]'F2 >= target_TIF3[1]-.1)\n",
"@constraint(mod2, target_TIF3[2]+.1 >= [fill(0,60);IIF3[2,:];fill(0,60)]'F2 >= target_TIF3[2]-.1)\n",
"@constraint(mod2, target_TIF3[3]+.1 >= [fill(0,60);IIF3[3,:];fill(0,60)]'F2 >= target_TIF3[3]-.1)\n",
"@constraint(mod2, target_TIF3[4]+.1 >= [fill(0,60);IIF3[4,:];fill(0,60)]'F2 >= target_TIF3[4]-.1)\n",
"@constraint(mod2, target_TIF3[5]+.1 >= [fill(0,60);IIF3[5,:];fill(0,60)]'F2 >= target_TIF3[5]-.1)\n",
"\n",
"@constraint(mod2, target_TIF3[1]+.1 >= [fill(0,120);IIF3[1,:]]'F3 >= target_TIF3[1]-.1)\n",
"@constraint(mod2, target_TIF3[2]+.1 >= [fill(0,120);IIF3[2,:]]'F3 >= target_TIF3[2]-.1)\n",
"@constraint(mod2, target_TIF3[4]+.1 >= [fill(0,120);IIF3[4,:]]'F3 >= target_TIF3[4]-.1)\n",
"@constraint(mod2, target_TIF3[3]+.1 >= [fill(0,120);IIF3[3,:]]'F3 >= target_TIF3[3]-.1)\n",
"@constraint(mod2, target_TIF3[5]+.1 >= [fill(0,120);IIF3[5,:]]'F3 >= target_TIF3[5]-.1)\n",
"\n",
"# between test forms\n",
"for j in 1:60\n",
" @constraint(mod2, 1 >= [F1[j] + F2[j+60] + F3[j+120]][1] >= 0)\n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"最後のloopは冊子間に重複項目を許さない制約についてです。loopにしてやらないと`Range constraint is not supported for GenericAffExpr{Float64,VariableRef}`で怒られます。\n",
"\n",
"それでは，最適化を実行します。"
]
},
{
"cell_type": "code",
"execution_count": 79,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 27.991487 seconds (263.30 k allocations: 34.654 MiB, 0.05% gc time)\n"
]
}
],
"source": [
"@time optimize!(mod2)"
]
},
{
"cell_type": "code",
"execution_count": 80,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"OPTIMAL::TerminationStatusCode = 1"
]
},
"execution_count": 80,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"termination_status(mod2) # OPTIMAL::TerminationStatusCode = 1 ならいちおうOK"
]
},
{
"cell_type": "code",
"execution_count": 91,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"3×10 Array{Int64,2}:\n",
" 1 8 29 32 33 34 35 41 55 60\n",
" 3 10 11 12 26 28 30 37 48 59\n",
" 4 7 9 16 19 20 25 36 47 56"
]
},
"execution_count": 91,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"get_F1 = value.(F1)\n",
"get_F2 = value.(F2)\n",
"get_F3 = value.(F3)\n",
"\n",
"[([1:60; 1:60; 1:60][isone.(get_F1)])';([1:60; 1:60; 1:60][isone.(get_F2)])'; ([1:60; 1:60; 1:60][isone.(get_F3)])']"
]
},
{
"cell_type": "code",
"execution_count": 95,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 95,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot(tif([-4:0.1:4;], [α2 α2 α2][isone.(get_F1)], [β2 β2 β2][isone.(get_F1)]), label = \"optimized test F1\")\n",
"plot!(tif([-4:0.1:4;], [α2 α2 α2][isone.(get_F2)], [β2 β2 β2][isone.(get_F2)]), label = \"optimized test F2\")\n",
"plot!(tif([-4:0.1:4;], [α2 α2 α2][isone.(get_F3)], [β2 β2 β2][isone.(get_F3)]), label = \"optimized test F3\")\n",
"plot!(tif([-4:0.1:4;], α2, β2) * 10/60, label = \"target TIF\") # SUCCESS...?"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Julia 1.1.0",
"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
}