非線形計画法　練習プログラム　XNLP ver.1.0.2



０．目次

　１．はじめに
　２．ダウンロード＆インストール
　３．実行
　４．目的関数の書き方
　５．初期値の与え方
　６．エラーと実行結果
　７．最適化の方法と設定
　８．簡単な問題のリスト
　９．制約条件の付け方
１０．Ramseyモデルとその結果

１．はじめに

ＧＲＡＳで用いている最適化エンジンの本体を使った、非線形計画法の簡単な練習プログラムです。
簡単なといっても、１０期や２０期のRamseyモデルなど楽に解く力を持っているはずです。
変数１００個くらいは問題ないでしょう。

ＧＲＡSと決定的に違うことは、モデルの構造を自分で組むことができることです。ＧＲＡSはモデルの構造は
すでに組み込まれていて、シミュレーションはパラメータを変えるだけですが、このxnlpは、ちょっと面倒ですが、
自分で組むことができます。

変数を含んだ数式で与えて、それをxnlpが解釈して実行します。これが、私が前から組みたかったプログラムです。

まず、ダウンロードしてインストールしてください。それからまた解説します。

２．ダウンロード＆インストール

（１）下記のセットアップファイルをデスクトップなりに保存して、実行させてください。

XnlpSetup102 (2008年2月18日更新)
（無料ですが、無保証ですので自己責任で使ってください。）

マイクロソフトの.NET Framework ver.2（無料） がインストールされていない場合、
まず.NET Framework ver.2をインストールすることが要求されます。（すでにＧＲＡＳがインストールされていれば、
.NET Framework ver.2はインストール済みのはずですから、簡単に終わります。）

セットアップの途中で、Microsoft .NET FrameWork 2 が必要ですと、言われた場合には、
マイクロソフトのサイトからダウンロードしてインストールしてください。

.NET Framework ver.2よりも新しいバージョンver.3.0（Windows VISTAにインストール済み） ver.35でもかまいません。

（２）アンインストール
Xnlp をアンインストールするのには、二つの方法があります。
A. セットアップファイルをもう一度起動して、最初の画面で「削除」を選択する。
B. Windowsのコントロールパネルのプログラムの追加と削除からXnlpを選択し「削除」を実行する。
通常これらの方法で問題なく削除できます。

（※セキュリティソフトなどの影響で、時間がかかるなどの、問題が発生する場合があります。自己責任で対応願います。 ）

３．実行

デスクトップにできたアイコンをクリックすると、xnlpはスタートします。３つのテキストエリアをもったWindowフレームが
開きます。上から、目的関数入力エリア、初期値入力エリア、コンソール出力です。



簡単な問題を解かせてみましょう。
目的関数エリアに、

(x-1)^2+(y-2)^2

と入力してください。（コピーアンドペーストでＯＫ）。この式で、「x, y」は変数、「^」はべき乗の記号です。

目的関数は、x=1, x=2で最小解0を達成するのは一目瞭然ですね。（なお、xnlpは最小化問題しか解きません。
最大化問題を解かせたいときは、目的関数にマイナスをつけて最小化問題に変換してください。）

次に、初期値を与えなければなりません。初期値欄に次の文字列を入力してください。

1,1

これは、目的関数に登場する変数の、登場順に初期値をカンマで区切って与えます。カンマ以外のもので区切っては、
いけません。

続いて、フレーム右上のRUNボタンを押せば、瞬間的に答えが出ます。しかし、人間の目で自明でも、コンピュータは
目と脳をもっていませんので、初期値から本の狭い範囲の状況を見ながら、最小値に近づいていく手続きをとります。
ちょうど、地図を持って大きな山を登るようなもので、頂上は簡単には見いだせないのです。

次に、いきなり複雑な問題を解かせてみましょう。6期間の動学的経済モデル（Ramseyモデル）です。目的関数欄に、
つぎの文字列をコピー＆ペーストしてください。

-(c1^(1-0.7)/(1-0.7))-(1/1.1)*(c2^(1-0.7)/(1-0.7))
-(1/1.1)^2*(c3^(1-0.7)/(1-0.7))-(1/1.1)^3*(c4^(1-0.7)/(1-0.7))
-(1/1.1)^4*(c5^(1-0.7)/(1-0.7))-(1/1.1)^5*(c6^(1-0.7)/(1-0.7))
+10000*(1.3*10^0.6-c1-i1)^2
+10000*(1.3*k2^0.6-c2-i2)^2
+10000*(1.3*k3^0.6-c3-i3)^2
+10000*(1.3*k4^0.6-c4-i4)^2
+10000*(1.3*k5^0.6-c5-i5)^2
+10000*(1.3*k6^0.6-c6-i6)^2
+10000*(k2-i1-10)^2
+10000*(k3-i2-k2)^2
+10000*(k4-i3-k3)^2
+10000*(k5-i4-k4)^2
+10000*(k6-i5-k5)^2
+10000*(30-i6-k6)^2

ちょっと長くて意味不明かもしれませんが、のちに解説します。10000という数字は、制約条件を満たさないことに
対するペナルティです。

続いて、初期値の欄に次のように入力してください。

1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1

先と同じようにRUNボタンをクリックすると、私のパソコンの場合、瞬く間に答えを出します。
コンソールに初期設定などの値の表示が終わると次のような答えが出力されます。

Optimal Solution!
iteration = 213
ObjectValue = -23.47526427
c1 = 2.34391066
c2 = 2.91418448
c3 = 3.52654984
c4 = 4.17803129
c5 = 4.86855381
c6 = 5.58595450
i1 = 2.83151010
k2 = 12.83153764
i2 = 3.09633663
k3 = 15.92789577
i3 = 3.31633895
k4 = 19.24425182
i4 = 3.48719100
k5 = 22.73145663
i5 = 3.60216548
k6 = 26.33363338
i6 = 3.66635730

このあとに、
Encountered following errors, but it may not affect the result.
Error : 6 : Power Arith. Error
というメッセージが出ますが、この場合は、無視してかまいません。

c1～c6は第１期から第６期までの最適消費水準です。i1～i6は同じく投資水準。k2～k6は同じく投資水準。
k1は初期値として、10という値を与えていることが目的関数から読み取れるでしょう。第７期に30単位の
資本ストックを残さなければなりません。

では、この解が 本当に最適かといわれると困るのですが、ＧＡＭSという著名なソフトで同じ問題を解くと、
ほぼ同じ答えがでますので、間違いはないでしょう。小数点３ケタまでは一致します。もっと精密な解を
求めたければペナルティの値を大きくすればよいでしょう。

というわけで、一通りの実行練習は終わることにします。

４．目的関数の書き方

目的関数の書き方を説明します。

（１）目的関数には、すべて半角文字で、全文字列が10000字以内で記入しなければなりません。
たとえば、１行80文字で書いた場合は、125行までしか書けないということになります。

（２）目的関数内では、空白、タブ、改行などは読み飛ばされます。項の区切りとしての意味はあります。
すなわち、xvarという変数が、xv arと書かれたら別物として扱われます。したがって、時々
改行を入れて関数を見やすくすることができます。

（３）変数は、２５５文字以内でアルファベットで始まる、アルファベットか数字でできた文字列です。
新たな変数が登場する順で、登録していきます。

（４）変数の数の上限は、1000個にしてありますが、これほど使うことはまずあり得ないでしょう。

（５）演算子としては、+と-は普通のとおりです。掛け算は、*を使います。割り算は、/です。
べき乗は、＾を使います。

（６）カッコ( )は、学校で習った通りの使い方ができます。すなわち、その括弧の中の式が最優先で
評価されます。

（７）数学関数としては、次のものが使えます。
　　"sin", "cos", "tan", "sqrt", "log", "exp","abs"
　　左から、サイン、コサイン、タンジェント、平方根、自然対数、自然対数の底のべき乗、絶対値です。
　原則的に、それぞれの値あるいは変数は直後にカッコでくくって与えなければなりません。たとえば、
　log(x) といった調子です。しかし、log x でもxが単項として自立していれば、適切に処理しますが、
　これは裏技であり間違いのもとです。基本的に前者、log(x)という書き方に徹してください。

（８）演算の優先順位は、数値や数学関数が最優先、つづいてべき乗、さらに続いて、掛け算と割り算、
最後に、加減演算です。次の項につくマイナス記号も最優先項目になります。すなわち、-2^2はプラスの
4であり、マイナス4ではありません。演算順序がわからなくなったら、( )を使って、明示してください。
たとえば、-2^2は、(-2)^2と書けば、自分で確信を持つことができるでしょう。

５．初期値の与え方

初期値の与え方は、目的関数に登場する変数について、それが出てきた順で与える必要があります。
初期値は、半角のカンマで区切られている必要があります。それらしい初期値を与えたほうが、
xnlpは気持ちよく答えを出す可能性があります。

順序はチェックしませんが、与えた初期値の数と変数の数はぴったり一致する必要があります。
あわなければ、適当な値を入れておくことも可能だったのですが、エラーとしたほうが、
目的関数の間違いに気付きやすいので、厳密にチェックするようにしています。

６．エラーと実行結果

（１）目的関数のコンパイルと実行
目的関数がそのように書かれていなかったら、そのようにエラーを吐き出します。xnlpは、目的関数を
一度、コンパイルして中間言語に書き直して使っています。そのコンパイル段階のエラーは文法に
かかわるもので、実行段階のエラーは、負の値に対してルートをとるなど計算エラーになります。

（２）計算エラー
計算結果の後に
Encountered following errors, but it may not affect the result.
Error : 6 : Power Arith. Error
などと出力されることがあります。これは、計算途中でエラーがあったのですがxnlpがわで、適当に
エラーを回避したことを表します。このエラーがあっても、最適解を求められる場合が多いです。
しかし、まず、しっかり解をチェックしてください。また、このエラーを回避する処置をしたために、
解が非合理的になってしまう場合もあります。その場合は、目的関数を改良する必要が
あります。

（３）エラーメッセージ

コンパイル時と計算時に出るxnlpのエラーメッセージは次の通りです。

７．最適化の方法と設定

うまく収束しない場合も多々ありますので、その場合には、方法や設定を変えていろいろ試みる必要があります。

（１）基本的な最適化法
　xnlpの基本的な最適化法は、準ニュートン法です。GRASのサイトに若干の説明があるかもしれません。その過程で必要な
１次元探索（１次元最適化、直線探索とも言う）はArmijoの方法と、Wolfeの方法が使えます。トップフレーム上部のボタンで
切り替えることができます。一方で収束しない場合は、別の方法に切り替えてみてください。

（２）[Config]ボタンをクリックすると、詳しい設定ができるようになっています。これらについての説明もGRASにありますので
そちらを参考にしてください。

８．簡単な問題のリスト

簡単な問題のリストをいかに掲げています。ためしにやってみてください。必ずしも収束しません。局所解に陥ったり
するかもしれません。

(sin(x)-0.5)^2+(cos(y)-0.5)^2

((sin(x)-0.5)^2+(cos(y)-0.5)^2)+y/x
(１次元最適化法をWolfeに変更)

((sin(x)/y-0.5)^2+(cos(y)/x-0.5)^2)

(sqrt(sin(x))-0.5)^2+(cos(y)-0.5)^2

(x-1)^2+(y-2)^2

(sqrt(x)-0.5)^2+(sqrt(y)-0.8)^2

(x-0.5)^0.5+(y-0.8)^0.5+100*(xp^2-x+0.5)^2+100*(yp^2-y+0.5)^2
（なぜかうまく解けない、ちょっと無理があるか？）

(abs(x-0.5))^0.5+(abs(y-0.8))^0.5
(１次元最適化法をWolfeに変更)

-abs(x)^1*abs(y)^1+10*abs(x+2*y-10)
(１次元最適化法をWolfeに変更)

-0.2*log(x)-0.8*log(y)+1000000*(x+2*y-10)^2
（制約のもとでの効用最大化）

９．制約条件のつけ方

最初に示したラムゼーモデルの例を見ていただければ、制約条件の付け方は何となく、わかったと思います。

（１）等式制約の場合
たとえば、本来の目的関数に対して、2.0x+3.2y=10.0という制約条件を付けたければ、目的関数に加えて

1000*(2.0x+3.2y-10.0)^2

というペナルティ項を付け加えればいいわけです。ペナルティの付け方については、ほかにもいろいろあります。
1000という数字は、ペナルティのウェイトです。大きくなるほど、制約条件の厳密性を強く求めることになります。
あんまり大きくすると、本来の目的関数とのバランスが取れなくなります。

２乗にしなくて、たとえば、絶対値関数を使って

1000*abs(2.0x+3.2y-10.0)

でもよいようですが、理論的条件の関係で結果は必ずしもよいものにはなりません。

（２）不等式制約の場合
たとえば、2.0x+3.2y <= 10.0という制約条件をつけたい場合。ここで、"<="は"<"か"="のいずれかを満たすという
意味。このときは、新たな変数をスラック変数として、導入します。たとえばスラック変数をsとしましょう。
そして、この式を次のように等式条件に変換します。

2.0x+3.2y+s^2=10.0

この式は、かならず s^2 >= 0 であるから、上記の不等式条件と等価です。そして、先の等式条件と同じように

1000*(2.0x+3.2y+s^2-10.0)^2

として、目的関数に加えればよいわけです。不等式条件の、不等号が逆の場合どうすればよいかは、
練習問題とします。

何か、とても場当たり的な対応ですが、GRASで行っている高度な制約条件処理も、このような方法を
より理論的整合性を持たせる形で、かつ合理的にウェイトを調整するように行っているだけです。


１０．Ramseyモデルとその解き方

冒頭で示した例である、Ramseyモデルを少しだけ解説しよう。この型のモデルの３期バージョンについての
別の説明が、ブログ温暖化対策の経済学に説明されている。参考にしていただきたい。以下の例は、
それを６期間に伸ばしているだけである。

もう一度、その目的関数を示しておく。

-(c1^(1-0.7)/(1-0.7))-(1/1.1)*(c2^(1-0.7)/(1-0.7))
-(1/1.1)^2*(c3^(1-0.7)/(1-0.7))-(1/1.1)^3*(c4^(1-0.7)/(1-0.7))
-(1/1.1)^4*(c5^(1-0.7)/(1-0.7))-(1/1.1)^5*(c6^(1-0.7)/(1-0.7))
+10000*(1.3*10^0.6-c1-i1)^2
+10000*(1.3*k2^0.6-c2-i2)^2
+10000*(1.3*k3^0.6-c3-i3)^2
+10000*(1.3*k4^0.6-c4-i4)^2
+10000*(1.3*k5^0.6-c5-i5)^2
+10000*(1.3*k6^0.6-c6-i6)^2
+10000*(k2-i1-10)^2
+10000*(k3-i2-k2)^2
+10000*(k4-i3-k3)^2
+10000*(k5-i4-k4)^2
+10000*(k6-i5-k5)^2
+10000*(30-i6-k6)^2

まず、最初の３行はオリジナルの目的関数である、６期分の消費の流列によって与えられる総効用を
最大にしている。効用の割引率は10％に設定している。1/1.1=1/(1+0.1)は割引ファクターである。１期あたりの効用は、
限界効用の弾力性が0.7で一定しているものを用いている。

次の６行は、各期の需給バランス式である。すなわち、t期の需給バランスが

1.3*k(t)^0.6-c(t)-i(t) t=1,2,3,4,5,6

としてあらわされている。1.3は全要素生産性、0.6は資本の生産に対する弾力性であり、生産関数を
表している。

最後の６行は資本蓄積方程式である。

割引率を変えたり、いろいろ試みられたい。


以上