システムとモデリング

modelica, Julia, Design Structure Matrix, SysML, 他モデリング全般について。

【Julia言語おすすめパッケージ】Juliaの数式処理システム[Symbolics.jl]で数学の問題を解いてみる

今回も前回に引き続き数式処理システムSymbolics.jlの解説です。今回はこのパッケージを用いて簡単な多変数最適化の文章題を解いてみたいと思います。

問題は以下書籍からそのまま抜き出しました。

また、前回の記事は以下です。

otepipi.hatenablog.com

問題内容

2種類のテレビの最適生産量を求める問題です。

あるテレビメーカーが2つの新製品を計画している。ひとつは19インチLCDフラットパネルの製品で、メーカー希望小売価格は339ドルである。もう一つは21インチLCDフラットパネルの製品で、メーカー希望小売価格は399ドルである。メーカーにおけるコストは19インチの製品で1台195ドル、21インチの製品で1台225ドル、そしてその他に固定コストが400,000ドルである。これらの製品が販売される競争市場では年間売上台数が平均販売価格に影響を与える。いずれの製品も、1台多く売れるごとに平均販売価格が1セント下がると見積もられている。さらに、19インチの製品の販売は21インチの製品の販売に影響を与え、また逆の影響もある。それは、21インチの製品が1台多く売れるごとに19インチの製品の販売価格が0.3セント下がり、19インチの製品が1台多く売れるごとに、21インチの製品の平均販売価格が0.4セント下がると見積もられている。この製品はそれぞれ何台生産すべきだろうか?

利益を最大にする21インチ、19インチそれぞれの生産量を求めることになります。 では、この問題を解いていきたいと思います。

変数

問題文の中から読み取れる変数を整理します。

記号 意味
s 19インチの製品の販売台数(年あたり)
t 21インチの製品の販売台数(年あたり)
p 19インチの製品の販売価格(ドル)
q 21インチの製品の販売価格(ドル)
C 製造コスト(ドル/年)
R 製品販売による収入(ドル/年)
P 製品販売による利益(ドル/年)

Pを最大にするs,tの値はいくらか?が問題の趣旨となります。

仮定

問題文の中から読み取れる各記号の関係式は以下になります。

 p=339 - 0.01s - 0.003t

 q=399 - 0.004s - 0.01t

 R = ps + qt

 C = 400,000 + 195s + 225t

 P = R - C

 s ≧ 0

 t ≧ 0

Juliaで問題を解く

では、解いていきましょう。

パッケージの宣言

まずは何はさておきSymbolics.jlを宣言します。

using Symbolics

変数の宣言と19インチの販売価格の立式

変数s,tを宣言し、それらを用いて19インチの販売価格pとの関係式を宣言します。

@variables s t
p = 339 - 0.01s - 0.003t

\begin{equation} 339 - 0.01 s - 0.003 t \end{equation}

21インチの販売価格の立式

同様にs,tを用いて21インチの販売価格qを立式します。

q = 399 - 0.004s - 0.01t

\begin{equation} 399 - 0.004 s - 0.01 t \end{equation}

売上の計算

収入Rを求めます。販売価格×生産数量ですね。

R = p*s + q*t

\begin{equation} s \left( 339 - 0.01 s - 0.003 t \right) + t \left( 399 - 0.004 s - 0.01 t \right) \end{equation}

コストの計算

コストCを求めます。

C = 400000 + 195s + 225t

\begin{equation} 400000 + 195 s + 225 t \end{equation}

利益の計算

最大化したい利益Pを立式します。

P = R - C

\begin{equation} -400000 - 195 s + s \left( 339 - 0.01 s - 0.003 t \right) - 225 t + t \left( 399 - 0.004 s - 0.01 t \right) \end{equation}

展開してみます。

func = simplify(P;expand=true)

\begin{equation} -400000 + 144 s + 174 t - 0.01 s^{2} - 0.01 t^{2} - 0.007 s t \end{equation}

上記を最大にするs,tの値を求めます。 最大化したい関数を f(s,t)とすると

 ∇f = 0

となるs,tがこの問題の解になります。すなわち

 \frac{\partial f}{\partial s} = 0

 \frac{\partial f}{\partial t} = 0

上記を満たすs,tが解となります。これを求めてみましょう。

 \frac{\partial f}{\partial s}を求める

以下でsの偏微分を定義します

Ds = Differential(s)
(::Differential) (generic function with 2 methods)
Ds(func)

\begin{equation} \mathrm{\frac{d}{d s}}\left( -400000 + 144 s + 174 t - 0.01 s^{2} - 0.01 t^{2} - 0.007 s t \right) \end{equation}

式を展開して微分されていることを確認します。

fs = expand_derivatives(Ds(func))

\begin{equation} 144 - 0.02 s - 0.007 t \end{equation}

 \frac{\partial f}{\partial t}を求める

同様の方法でtの偏微分を求めます。

Dt = Differential(t)
Dt(func)
ft = expand_derivatives(Dt(func))

\begin{equation} 174 - 0.007 s - 0.02 t \end{equation}

連立方程式を解く。

仕上げに下記連立方程式を解きます。

 \frac{\partial f}{\partial s} = 0

 \frac{\partial f}{\partial t} = 0

eqs = [
        fs ~ 0,
        ft ~ 0]
var = [s,t]
Symbolics.solve_for(eqs,var)

\begin{equation} \left[ \begin{array}{c} 4735.042735042735 \ 7042.735042735043 \ \end{array} \right] \end{equation}

結果、以下のように正答を求めることができました。

 19インチの生産量 = 4735台

 21インチの生産量 = 7043台

簡単ですが今回はここまでにします。