はじめに
Juliaは簡単に分数を扱えると聞いていて、実際にやってみたら確かに強かった、という記事です。
本記事では線形計画問題の基底を変更する操作をします。
扱う問題
$$
\begin{align*}
\text{minimize} ~ & ~ z = -x_1 -x_2\\
\text{subject to} ~ & ~ 5x_1 + 8x_2 + x_3 = 96\\
& ~ 7x_1 + 4x_2 + x_4 = 84\\
& ~ x_1, x_2, x_3, x_4 \geq 0
\end{align*}
$$
この問題の基底を$x_1$と$x_2$にします。
計算をするために式を行列にしておきます。
$$
x = \left(
\begin{array}{c}
x_1\\
x_2\\
x_3\\
x_4\\
\end{array}
\right), ~~
A = \left(
\begin{array}{cccc}
5 & 8 & 1 & 0\\
7 & 4 & 0 & 1\\
\end{array}
\right), ~~
b = \left(
\begin{array}{c}
96\\
84\\
\end{array}
\right), ~~
c = \left(
\begin{array}{c}
-1\\
-1\\
0\\
0\\
\end{array}
\right)\\
\begin{align*}
\text{minimize} ~ & ~ c^Tx\\
\text{subject to} ~ & ~ Ax = b\\
& ~ x \geq \bf{0}
\end{align*}
$$
$x_1, x_2$を基底にするので$A$の1,2列を取り出して$B$を定義します。
$$
B = \left(
\begin{array}{cc}
5 & 8\\
7 & 4\\
\end{array}
\right)
$$
Juliaでの計算
まずは計算に必要な行列などを一通り定義します。
b = [
96
84
]
A = [
5 8 1 0
7 4 0 1
]
c = [
-1
-1
0
0
]
B = [
5//1 8//1
7//1 4//1
]
Bの型はMatrix{Rational{Int64}}になります。
Bに関する計算では浮動小数ではなく分数で計算してくれます。
変形後の成約の部分$B^{-1}b$と$B^{-1}A$を計算します。
display(B \ b)
display(B \ A)
2-element Vector{Rational{Int64}}:
8//1
7//1
2×4 Matrix{Rational{Int64}}:
1//1 0//1 -1//9 2//9
0//1 1//1 7//36 -5//36
目的関数から基底変数を消去して、定数部分も計算します。
$c = (-1 ~ -1 ~~ 0 ~~ 0)^\text{T}$なので$B^{-1}A$の1行目と2行目に1をかけて足します。
display(c + vec(sum(B \ A, dims=1)))
display(-sum(B \ b))
4-element Vector{Rational{Int64}}:
0//1
0//1
1//12
1//12
-15//1
(inv(B)A)[1, 1:4]の部分は行列の1行目を取り出したものですが、これはベクトルVector{Rational{Int64}}型になります。行で取り出しても縦に要素が並ぶベクトルになるので少し注意が必要です。
コメントにて推奨いただいたコードに変更して、1行4列のMatrixを返すsumをvecでベクトルにするコードに変更しました。
変形後の形
Juliaでの計算で変形後の線形計画問題が書けます。
$$
\begin{align*}
\text{minimize} ~ & ~ \frac{1}{12}x_3 + \frac{1}{12}x_4 - 15\\
\text{subject to} ~ & ~ x_1 - \frac{1}{9}x_3 + \frac{2}{9}x_4 = 8\\
& ~ x_2 + \frac{7}{36}x_3 - \frac{5}{36}x_4 = 7\\
& ~ x_1, x_2, x_3, x_4 \geq 0
\end{align*}
$$
まとめ
B = [
5//1 8//1
7//1 4//1
]
のようにしてMatrix{Rational{Int64}}型の変数が定義できる。
バージョン等
- Docker Image jupyter/datascience-notebook 2eac078be835
Jupyter Labを使用 - Julia 1.6.0
参考ページ