optimize_print.book

356  Download (0)

Full text

(1)

For Use with MATLAB

®

Optimization

Toolbox

(2)

営業部 techmatlab@cybernet.co.jp 技術部 03-5978-5411 電話 03-5978-5440 ファクシミリ 東京都文京区大塚 住所 2 丁目 15 番地 6 号 ニッセイ音羽ビル

For contact information about worldwide offices, see the MathWorks Web site. Optimization Toolbox User’s Guide

 COPYRIGHT 1990 - 2004 by The MathWorks, Inc.

The software described in this document is furnished under a license agreement. The software may be used or copied only under the terms of the license agreement. No part of this manual may be photocopied or repro-duced in any form without prior written consent from The MathWorks, Inc.

FEDERAL ACQUISITION: This provision applies to all acquisitions of the Program and Documentation by, for, or through the federal government of the United States. By accepting delivery of the Program or Documentation, the government hereby agrees that this software or documentation qualifies as commercial computer software or commercial computer software documentation as such terms are used or defined in FAR 12.212, DFARS Part 227.72, and DFARS 252.227-7014. Accordingly, the terms and conditions of this Agreement and only those rights specified in this Agreement, shall pertain to and govern the use, modification, reproduction, release, performance, display, and disclosure of the Program and Documentation by the federal government (or other entity acquiring for or through the federal government) and shall supersede any conflicting contractual terms or conditions. If this License fails to meet the government's needs or is inconsistent in any respect with federal procurement law, the government agrees to return the Program and Documentation, unused, to The MathWorks, Inc.

MATLAB, Simulink, Stateflow, Handle Graphics, and Real-Time Workshop are registered trademarks, and TargetBox is a trademark of The MathWorks, Inc.

Other product or brand names are trademarks or registered trademarks of their respective holders.

Printing History: November 1990 First printing

December 1996 Second printing For MATLAB 5

January 1999 Third printing For Version 2 (Release 11) September 2000 Fourth printing For Version 2.1 (Release 12)

June 2001 Online only Revised for Version 2.1.1 (Release 12.1) September 2003 Online only Revised for Version 2.3 (Release 13SP1) June 2004 Fifth printing Revised for Version 3.0 (Release 14) October 2004 翻訳

本書の内容の一部あるいは全部を無断で転載、複製、複写することを禁じます。 本書の内容は予告なく変更することがあります。

(3)

Acknowledgements

The MathWorks would like to acknowledge the following contributors to the Optimization Toolbox.

Thomas F. Coleman researched and contributed the large-scale algorithms for constrained and unconstrained minimization, nonlinear least squares and curve fitting, constrained linear least squares, quadratic programming, and nonlinear equations.

Dr. Coleman is Professor of Computer Science and Applied Mathematics at Cornell University. He is Director of the Cornell Theory Center and the Cornell Computational Finance Institute. Dr. Coleman is Chair of the SIAM Activity Group on Optimization, and a member of the Editorial Boards of Applied

Mathematics Letters, SIAM Journal of Scientific Computing, Computational Optimization and Applications, Communications on Applied Nonlinear Analysis, and Mathematical Modeling and Scientific Computing.

Dr. Coleman has published 4 books and over 70 technical papers in the areas of continuous optimization and computational methods and tools for

large-scale problems.

Yin Zhang researched and contributed the large-scale linear programming algorithm.

Dr. Zhang is Associate Professor of Computational and Applied Mathematics on the faculty of the Keck Center for Computational Biology at Rice University. He is on the Editorial Board of SIAM Journal on Optimization, and is Associate Editor of Journal of Optimization: Theory and Applications.

Dr. Zhang has published over 40 technical papers in the areas of interior-point methods for linear programming and computation mathematical

(4)
(5)

Contents

1

はじめに

Optimization Toolboxとは ? . . . 2 最適化の例題 . . . 3 問題 . . . 3 問題の設定 . . . 3 解を探す . . . 4 その他の例題 . . . 5

2

チュートリアル

イントロダクション . . . 3 ツールボックスで対象とされる 問題 . . . 3 最適化関数の使用 . . . 6 中規模および大規模アルゴリズム . . . 7 標準的なアルゴリズムを使用した例題 . . . 9 制約なし最小化の例題 . . . 10 非線形不等式制約付きの例題 . . . 11 範囲指定を伴った制約付きの例題 . . . 13 勾配指定を伴った制約付きの例題 . . . 14 勾配のチェック : 解析的に求めた勾配と数値的に求めた勾配 . . . 16 等式制約付き例題 . . . 16 最大化 . . . 18 ゼロより大きい制約 . . . 18 無名関数と入れ子関数によるグローバル変数の回避 . . . 19 解析ヤコビ行列による非線形方程式 . . . 22 有限差分ヤコビ行列による非線形方程式 . . . 25 複数の目的関数をもつ例題 . . . 26 大規模例題 . . . 39

(6)

勾配およびヘッセ行列を使った非線形最小化 . . . 49 勾配および ヘッセ行列にスパースパターンを使った非線形最小化 50 範囲による制約および範囲付きの前処理演算子を使った非線形最小化 52 等式制約を使った非線形最小化 . . . 56 密であるが構造化されたヘッセ行列および等式制約を伴った非線形最小化 57 範囲による制約をもつ 2 次型の最小化 . . . 61 ヘッセ行列が密であるが構造化されている 2 次型の最小化 . . . 63 範囲による制約を課した線形最小二乗法 . . . 68 等式および不等式を使った線形計画法 . . . 69 等式の中の密の列をもつ線形計画法 . . . 70 デフォルト オプション設定 . . . 73 デフォルト設定の変更 . . . 73 繰り返し結果の出力表示 . . . 76 出力見出し : 中規模アルゴリズム . . . 76 出力見出し : 大規模最適化アルゴリズム . . . 79 出力関数を繰り返しコールする . . . 81 例題で行うこと . . . 81 出力関数 . . . 82 例題のための M- ファイルの作成 . . . 83 例題の実行 . . . 84 M-ファイルの代わりに無名関数を使った最適化 . . . 88 標準的な問題と対処法 . . . 90 参考文献 . . . 93

(7)

3

標準的なアルゴリズム

最適化の概要 . . . 3 制約なしの最適化 . . . 4 準 Newton 法 . . . 5 ライン探索法 . . . 7 準 Newton 法の実行 . . . 9 ヘッセ行列の更新 . . . 9 ライン探索法 . . . 9 最小二乗最適化 . . . 16 Gauss-Newton 法 . . . 17 Levenberg-Marquardt 法 . . . 18 非線形最小二乗法の実行 . . . 19 連立非線形方程式の解法 . . . 22 Gauss-Newton法 . . . 22 Trust-Region Dogleg 法 . . . 22 非線形方程式の実行 . . . 24 制約付き最適化 . . . 25 逐次二次計画法 (SQP) . . . 26 QP サブ問題 . . . 26 SQP 法の実行 . . . 28 シンプレックスアルゴリズム . . . 33 複素目的関数の最適化 . . . 39 多目的最適化の紹介 . . . 39 ゴール到達法 . . . 45 ゴール到達法のアルゴリズム実行 . . . 46 参考文献 . . . 49

(8)

非線形最小化に対する Trust-Region 法 . . . 2 前処理を用いた共役勾配法 . . . 5 線形制約付き問題 . . . 7 線形等号制約 . . . 7 ボックス制約 . . . 7 非線形最小二乗法 . . . 10 2次計画法 . . . 11 線形最小二乗 . . . 12 大規模線形計画法 . . . 13 メインアルゴリズム . . . 13 前処理 . . . 16 参考文献 . . . 17

5

関数リファレンス

関数カテゴリ順リスト . . . 2 最小化 . . . 3 方程式の解法 . . . 3 最小二乗法 ( カーブフィッティング ) . . . 3 ユーティリティ . . . 4 大規模法を使ったデモンストレーション . . . 4 中規模法を使ったデモンストレーション . . . 4 関数の引数 . . . 5 入力引数 . . . 5 出力引数 . . . 7

(9)

最適化オプション . . . 9 出力関数 . . . 14 関数− アルファベット順のリスト . . . 23

(10)
(11)

1

はじめに

Optimization Toolboxとは ? (p. 1-2) ツールボックスを紹介し、解くことができる問題のタイプ を述べます。

(12)

Optimization Toolbox

とは ?

Optimization Toolbox は、 MATLAB®の数値計算環境機能を拡張した関数を集めた

ものです。このツールボックスは、つぎに示す多くの最適化タイプのルーチンを 備えています。 • 制約なし非線形最小化 • ゴール到達問題、ミニマックス問題、および半無限最小化問題を含む制約付き 非線形最小化 • 二次および線形計画法 • 非線形最小二乗およびカーブフィッテリング • 非線形方程式システムの解法 • 制約付き線形最小二乗法 • スパースおよび構造化された大規模問題 ツールボックス関数は、すべて、特定の最適化アルゴリズムを実行する MATLAB ステートメントにより作成された MATLAB の M- ファイルです。つぎのステート メントにより、これらの関数の MATLAB コードを見ることができます。 type function_name ユーザの M- ファイルを記述したり、他のツールボックスや MATLAB、Simulink® などと組み合わせて使用することにより、Optimization Toolbox の機能を拡張する ことができます。

(13)

最適化の例題

最適化の例題

この節では、ツールボックスの線形最小二乗問題を解く関数 lsqlin を使用して 最適化問題を解く方法を説明する例題を述べます。この節では、つぎのトピック スを述べます。 • 1-3ページの " 問題 " • 1-3ページの " 問題の設定 " • 1-4ページの " 解を探す " • 1-5ページの " その他の例題 "

問題

この例題の問題は、平面 上の原点に最も近い点を探すことで す。この問題を解くための最も簡単な方法は、平面上の点 から 原点への距離の二乗を最小化することです。この場合も、実際の距離を最小化す るのと同じ最適化点が出力されます。任意の点 から原点への距離の 二乗は、 なので、この問題をつぎのように記述することができます。 制約条件 関数 f(x) は目的関数と呼ばれ、 は、等式制約です。より複雑 な問題では、他の等式制約、不等式制約、上限または下限の制約を含むことがあ ります。

問題の設定

この節では、つぎの形式の線形最小二乗問題を解く関数 lsqlinを適用する前に、 問題を設定する方法を示します。 x1+2x2+4x3 = 7 x = (x1, ,x2 x3) x1, ,x2 x3 ( ) x12+x22+x32 minimize x f x( ) x1 2 x22 x32 + + = x1+2x2+4x3 = 7 x1+2x2+4x3 = 7 minimize x f x( ) Cxd 2 =

(14)

ここで、 は、 Cx - d のノルムの二乗で、制約

があります。問題を設定するためには、パラメータ C, d, A, b, Aeq, beq に対して変

数を作成する必要があります。lsqlinは、つぎのシンタックスを使用して、入力

引数としてこれらの変数を受け取ります。 x = lsqlin(C, d, A, b, Aeq, beq)

変数を作成するためには、つぎのステップを行います。

1.

目的関数に対する変数の作成

を最小化したいため、Cx - d = x となるように C を 3 × 3 単 位行列、 d を 3 × 1 ゼロベクトルに設定できます。 C = eye(3); d = zeros(3,1);

2.

制約に対する変数の作成

この例題は、不等式制約をもたないため、入力引数において Aおよび bを空の行 列に設定できます。 等式制約 を、行列形式で

と表すことができます。ここで、 Aeq = [1 2 4] beq = [7]です。 Aeq および beq に対

する変数を作成するためには、つぎのように入力します。 Aeq = [1 2 4];

beq = [7];

解を探す

最適化問題を解くためには、

[x, fval] =lsqlin(C, d, [], [], Aeq, beq)

Cxd 2 Axb Aeq x⋅ = beq x12+x22+x32 = x 2 x1+2x2+4x3 = 7 Aeq x⋅ = beq

(15)

最適化の例題 と入力します。lsqlinは、つぎのように出力します。 x = 0.3333 0.6667 1.3333 fval = 2.3333 最小は点 xで最小となり、fvalは xから原点への距離の二乗です。 注意 この例題では、lsqlinは、デフォルトの 大規模アルゴリズム から 中規模 アルゴリズムに切り替えているという警告を出します。このメッセージは、結果 に何も影響しないため、無視しても問題はありません。 2-6ページの " 最適化関数 の使用 " は、大規模および中規模アルゴリズムについて、より詳細を提供します。

その他の例題

つぎの節では、最適化問題を解くその他の例題を述べます。 • 2-9ページの " 標準的なアルゴリズムを使用した例題 " • 2-39ページの " 大規模例題 "

(16)
(17)

2

チュートリアル

チュートリアルの章では、ツールボックス関数の使用方法についての情報を提供します。また、種々の 最適化問題を解くための例題を提供します。次の節からなります。 イントロダクション (p. 2-3) 最小化、方程式解法および最小二乗法、またはデータ フィッティング問題を解くことのできる関数を表形式で要 約します。また、中規模および大規模問題で使用可能なア ルゴリズムとライン探索法についての紹介と最適化ルーチ ンの使用に対する基本的なガイドラインを与えます。 標準的なアルゴリズムを使用した例題 (p. 2-9) 最小化の例題の節では、中規模アルゴリズムを示します。 これらの例では、ユーザ設定の勾配のない、またはある問 題と同様に、制約のない、または制約付きの問題を含みま す。また、この節では、最大化、0 より大きい制約、付加的 な引数、および複数の目的関数の例題も論じます。 大規模例題 (p. 2-39) 大規模問題例題の節では、大規模アルゴリズムを示しま す。これらの例では、制約のない、または制約付きの問題 を含む、スパース性の構造、前処理演算子の指定を含みま す。 デフォルト オプション設定 (p. 2-73) デフォルトパラメータ設定の使用と、変更方法について説 明します。また、関数を指定して、パラメータを定義する 方法を示し、一般に使用されるいくつかのパラメータの設 定例を提供します。 繰り返し結果の出力表示 (p. 2-76) 中規模、および大規模アルゴリズムの両方の繰り返し結果 で使用される列見出しについて示します。 出力関数を繰り返しコールする (p. 2-81) 各繰り返しで最適化関数に出力関数をコールさせる方法を 述べます。 M-ファイルの代わりに無名関数を使っ た最適化 (p. 2-88) 文字を含む式から無名関数を作成することにより、コマン ドラインで数学関数をあらわす方法を示します。

(18)

標準的な問題と対処法 (p. 2-90) アルゴリズムを効率的に改善し、共通の複雑な問題をうま くクリアし、そして標準的でない問題を変換することで、 使用する最適化関数の解法を改善するための助言を与えま す。 参考文献 (p. 2-93) Optimization Toolbox の概念をサポートしている文献のリス トです。

(19)

イントロダクション

イントロダクション

最適化とは、目的関数と通常呼ばれる関数の最小あるいは最大を見つける過程で す。 Optimization Toolbox は、一般的な非線形関数について最小化 ( あるいは最大 化 ) を実行する関数から構成されています。また、非線形方程式解法および最小 二乗 ( データフィッティング ) 問題のための関数も備えています。 このイントロダクションは、以下の節からなります。 • ツールボックスで対象とされる 問題 • 最適化関数の使用

ツールボックスで対象とされる 問題

つぎの表は、最小化、方程式解法、および最小二乗またはデータフィッティング 問題を解くために使用可能な関数を示したものです。 注意 下記の表は、問題が複雑になる順番にリストしています。 表 2-1: 最小化 タイプ 記法 関数 スカラー最小化 条件 fminbnd 制約なし最小化 fminunc,fmin search 線形計画法 条件 linprog 2次計画法 条件 quadprog f a( ) a min a1≤ ≤a a2 f x( ) x min fTx x min A x⋅ ≤b, Aeq x⋅ = beq, l≤ ≤x u 1 2 ---xTHx+fTx x min A x⋅ ≤b, Aeq x⋅ = beq, l≤ ≤x u

(20)

制約付き最小化 条件 fmincon ゴール到達 条件 fgoalattain ミニマックス 条件 fminimax 半無限最小化 条件 fseminf 表 2-1: 最小化 ( 続き ) タイプ 記法 関数 f x( ) x min c x( )≤0, ceq x( )= 0 A x⋅ ≤b, Aeq x⋅ =beq, l≤ ≤x u γ x,γ min F x( )–wγ ≤ goal c x( )≤0, ceq x( )= 0 A x⋅ ≤b, Aeq x⋅ =beq, l≤ ≤x u min x maxF i { } {Fi( )x } c x( )≤0, ceq x( )= 0 A x⋅ ≤b, Aeq x⋅ =beq, l≤ ≤x u f x( ) x min K x w( , )≤0 for all w c x( )≤0, ceq x( )= 0 A x⋅ ≤b, Aeq x⋅ =beq, l≤ ≤x u 表 2-2: 方程式の解法 タイプ 記法 関数 線形方程式 , n 個の方程式、n 個の変数 \ ( スラッ シュ ) 1変数による非線形方程式 fzero C x⋅ = d f a( ) = 0

(21)

イントロダクション 非線形方程式 , n 個の方程式、n 個の変数 fsolve 表 2-2: 方程式の解法 ( 続き ) タイプ 記法 関数 F x( ) = 0 表 2-3: 最小二乗法 ( カーブフィッティング ) タイプ 記法 関数 線形最小二乗法 , m 個の方程式 , n 個の変数 \ ( スラッ シュ ) 非負線形最小二乗法 条件 lsqnonneg 制約付き線形最小二乗法 条件 lsqlin 非線形最小二乗法 条件 lsqnonlin 非線形カーブフィッティ ング 条件 lsqcurvefit C x⋅ –d 2 x min 2 C x⋅ –d 2 x min 2 x≥0 C x⋅ –d 2 x min 2 A x⋅ ≤b, Aeq x⋅ =beq, l≤ ≤x u 1 2 --- F x( )22 x min 1 2 --- Fi( )x2 i

= l≤ ≤x u 1 2 --- F x xdata( , )–ydata 22 x min l≤ ≤x u

(22)

最適化関数の使用

この節では、最適化関数の使用についての基本的な情報を提供します。

目的関数の定義

最適化関数の多くでは、目的関数を計算する MATLAB 関数の作成が必要になり ます。この関数は、ベクトル入力を受け取り、タイプdoubleのスカラー出力を返 す必要があります。目的関数の作成方法は、2 通りあります。 • コマンドラインで無名関数を作成します。たとえば、 x2に対して無名関数を作 成するためには、つぎのように入力します。 square = @(x) x.^2; この後、第一入力引数としてsquareを使用して最適化関数を呼びます。この 方法は、目的関数が比較的簡単であり、後の MATLAB セッションで再使用す る必要がない場合に使用できます。 • この関数に対する M- ファイルをかきます。たとえば、 M-ファイルとして関数 x2をかくためには、 MATLAB エディタで新しいファイルを開き、つぎのコー ドを入力します。 function y = square(x) y = x.^2; この後、最初の入力引数として、@squareを使用して最適化関数をコールする ことができます。記号 @は、squareに対する関数ハンドルを作成します。目的 関数が複雑であったり、2 つ以上の MATLAB セッションで使用する予定があ る場合は、この方法を使用してください。

注意 Optimization Toolbox の関数は、タイプ doubleの入力のみを受け取るため、

目的関数は、タイプ doubleのスカラーの出力を返す必要があります。

最大化と最小化

ツールボックスの最適化関数は、目的関数を最小化します。関数 f を最大化する

ためには、 -f を最小化するように最適化関数を適用します。f が最大となる結果の

(23)

イントロダクション

オプションの変更

関数optimsetを使用して、入力引数としてoptions構造体を渡すことより、最 適化関数のデフォルトのオプションを変更できます。詳細は、 2-73 ページの " デ フォルト オプション設定 " を参照してください。

勾配を与える

最適化関数の多くは、最小を探すために目的関数の勾配を使用します。勾配を計 算し、options構造体を用いて最適化関数に渡す関数をかくことができます。 2-14 ページの " 勾配指定を伴った制約付きの例題 " では、このための例題を提供しま す。勾配関数を与えると、最適化関数の精度と速度が改善されます。ただし、目 的関数によっては、勾配関数を提供することができず、この場合、最適化関数は、

adaptive finite-difference methodを使用して勾配を計算します。

中規模および大規模アルゴリズム

本ガイド では、" 中規模 " アルゴリズムと " 大規模 " アルゴリズムを別々に取り 扱っています。中規模は、一般的な用語ではありません。本書では、大規模問題を 効果的に扱うために設計されている大規模アルゴリズムとの違いを示すために使 用しています。

中規模アルゴリズム

Optimization Toolbox のルーチンでは、アルゴリズムとライン探索法を選択するこ とができます。制約なし最小化の主要アルゴリズムは、 Nelder-Mead シンプレック ス探索法および BFGS 準 Newton 法です。制約付き最小化、ミニマックス、ゴール 到達、および半無限最適化は、逐次二次計画法 (SQP) をベースにしています。非 線形最小二乗問題では、Gauss-Newton 法および Levenberg-Marquardt 法を使用し

ます。非線形方程式解法でも、trust region dogleg アルゴリズムを使用します。

制約なし最小化および非線形最小二乗問題に対して、ライン探索法を選択します。 ライン探索法には、保護 3 次と 2 次の内挿および外挿法が使用されます。

大規模最適化アルゴリズム

線形計画法を除くすべての大規模アルゴリズムは、 trust-region 法を使います。範 囲による制約付き問題を解くには、reflective Newton 法を使用します。等式制約付 き 問 題 を 解 く に は、射 影 的 前 処 理 を 用 い た 共 役 勾 配 繰 り 返 し 法 (projective

preconditioned conjugate gradient iteration) を使用します。スパース行列を対象にし た繰り返しソルバ、あるいはスパースを直接処理するソルバを使用してカレント

ステップを決定し、線形システムを解くことができます。また、繰り返しソルバ

(24)

線形計画法は、 Mehrotra の予測子修正子アルゴリズムを変形した、主双対内点法 です。

(25)

標準的なアルゴリズムを使用した例題

標準的なアルゴリズムを使用した例題

本チュートリアルでは、Optimization Toolbox の中規模 ( つまり、標準的な ) アル

ゴリズムを示します。チュートリアルの最初の部分と似た例 (" 制約なし最小化の

例題 " から " 等式制約付き例題 " まで ) は、 M-ファイルoptdemoの 1 つ目のデモ

ンストレーション Tutorial Walk Through にも見られます。

注意 中規模という言葉は一般的な用語ではなく、 4-1ページの " 大規模アルゴリ

ズム " の章で説明されている大規模アルゴリズムと区別するために使用されてい る言葉です。

チュートリアルでは、関数 fminunc, fmincon, fsolveを使用します。その他の最

適化ルーチン、fgoalattain, fminimax, lsqnonlin, fseminfについてもほとん

ど同じ手法が使用されます。異なる点は問題の定式化と終了基準のみです。 2-26 ページの " 複数の目的関数をもつ例題 " では、多目的最適化を論じ、lsqnonlin, fminimax, fgoalattain を 使 用 す る い く つ か の 例 題 を 示 し ま す。こ れ に は、 Simulink をツールボックスと接続して使用する方法も含まれます。 この節には以下の例題が含まれます。 • 制約なし最小化の例題 • 非線形不等式制約付きの例題 • 範囲指定を伴った制約付きの例題 • 勾配指定を伴った制約付きの例題 • 勾配のチェック : 解析的に求めた勾配と数値的に求めた勾配 • 等式制約付き例題 以下についても論じます。 • 最大化 • ゼロより大きい制約 • 無名関数と入れ子関数によるグローバル変数の回避 • 解析ヤコビ行列による非線形方程式 • 有限差分ヤコビ行列による非線形方程式 • 複数の目的関数をもつ例題

(26)

制約なし最小化の例題

つぎの式の解の集合 [x1, x2] を求める問題を考えてみましょう。 (2-1) この 2 次元問題を解くために、関数値を戻す M- ファイルを記述します。そして、 制約なし最小化ルーチン fminuncを呼び出します。

ステップ 1: M- ファイル objfun.m の記述

function f = objfun(x) f = exp(x(1))*(4*x(1)^2+2*x(2)^2+4*x(1)*x(2)+2*x(2)+1);

ステップ 2: 制約なし最適化ルーチンの 1 つの呼び出し

x0 = [-1,1]; % 初期値推定 options = optimset('LargeScale','off'); 40回の関数計算の後、つぎの解を出します。 x = 0.5000 -1.0000 解xでの関数値がfvalに戻されます。 fval = 1.3030e-10 exitflagは、アルゴリズムが収束したかどうかを示します。exitflag > 0 は、極 小値が求まったことを意味します。 exitflag = 1 構造体outputには、最適化に関する詳細が与えられます。これは、fminuncに関

しては、iterations に繰り返し回数、 funcCountに関数値計算回数、stepsize

に最終のステップサイズ、firstorderoptに 1 次最適性の尺度 ( 制約なしの場合、 解における勾配の無限ノルムになります ) 、および algorithmに使用アルゴリズ ムタイプを含んだものになります。 output = iterations: 7 minimize x f x( ) e x1 4x21+2x22+4x1x2+2x2+1 ( ) =

(27)

標準的なアルゴリズムを使用した例題

funcCount: 40 stepsize: 1

firstorderopt: 9.2801e-004

algorithm: 'medium-scale: Quasi-Newton line search'

極小値が複数ある場合、ベクトル [x1,x2] に対する初期推定が、関数の計算回数と 解の値の両方に影響を与えます。上記の例題では、x0は[-1,1] に初期化されて います。 変数 optionsを fminuncに渡して、最適化アルゴリズムの特性を変更すること ができます。たとえば、つぎのようにします。 x = fminunc(@objfun,x0,options); optionsは、終了条件の許容範囲やアルゴリズムの選択に関する値を含む構造体 です。構造体 optionsは、関数optimsetを使って作成します。 options = optimset('LargeScale','off'); この例では、大規模アルゴリズムのデフォルト選択をやめ、中規模アルゴリズム を使用しました。他のオプションとしては、最適化繰り返し段階での途中経過の 表示頻度、終了条件に対する許容量、ユーザ設定の勾配またはヤコビ行列 のどち ら を 使 用 す る か、繰 り 返 し ま た は 関 数 値 計 算 の 最 大 回 数 等 が 含 ま れ ま す。 optimsetと個々の最適化関数については、 5-9ページの " 最適化オプション " を 参照してください。

非線形不等式制約付きの例題

不等式制約が(2-1)式 に付加されている場合、その結果の問題は、fminconを 使って解くことができます。たとえば、下記の制約のもとで、つぎの関数を最小 にする x を求めることを考えましょう。 (2-2) つぎの制約に従います。 いずれの制約も線形でないため、制約をコマンドラインでfminconに渡すことが できません。代わりに、2 つ目の M- ファイルconfun.mを作成して、ベクトル c の中に、カレント xにおける両方の制約の値を戻します。このとき、制約付き最 minimize x f x( ) e x1 4x21+2x22+4x1x2+2x2+1 ( ) = x1x2x1x2≤–1.5 x1x2≥–10

(28)

適子 fminconを呼び出します。fminconは、 の形式で記述される制約を 前提としているため、つぎのように制約を書き直す必要があります。

(2-3)

ステップ 1: 制約に対する M- ファイル confun.m の記述

function [c, ceq] = confun(x) % 非線形不等式制約 c = [1.5 + x(1)*x(2) - x(1) - x(2); -x(1)*x(2) - 10]; % 非線形等式制約 ceq=[];

ステップ 2: 制約付き最適化ルーチンの呼び出し

x0 = [-1,1]; % 解の初期推定値 options = optimset('LargeScale','off'); [x, fval] = ... fmincon(@objfun,x0,[],[],[],[],[],[],@confun,options) 38 回の関数値計算後、解 xを使って、関数値fvalは、つぎのようになります。 x = -9.5474 1.0474 fval = 0.0236 解 x の点での制約量も計算できます。 [c,ceq] = confun(x) これは、つぎの出力をします。 c= 1.0e-14 * 0.1110 -0.1776 ceq = [] c x( )≤0 x1x2x1x2+1.5≤0 x1x2–10 – ≤0

(29)

標準的なアルゴリズムを使用した例題 両制約値とも 0 以下であることに注意してください。すなわち、xは を満 たします。

範囲指定を伴った制約付きの例題

xの中の変数は、制約付き最適化関数に簡単な範囲による設定を行うことで、あ る範囲に限定することができます。fminconについて、コマンド x = fmincon(@objfun,x0,[],[],[],[],lb,ub,@confun,options); は、xを範囲lb <= x <= ubに制限します。 (2-2)式のx を 0 より大きいもの ( など ) に制限するには、つぎの コマンドを使用してください。 x0 = [-1,1]; % 解の初期推定値 lb = [0,0]; % 下限の設定 ub = [ ]; % 上限の設定なし options = optimset('LargeScale','off'); [x,fval = ... fmincon(@objfun,x0,[],[],[],[],lb,ub,@confun,options) [c, ceq] = confun(x) fminconへの 7 番目の引数として下限を設定するには、3 から 6 番目の値を指定 する必要があります。この例題では、線形不等式や線形等式制約をもっていない ので、これらの引数を []で指定しました。 13 回の関数計算の後、出力される解は、つぎの通りです。 x = 0 1.5000 fval = 8.5000 c= 0 -10 ceq = [] lbあるいは ubの要素が xより少ない場合、xの要素の内、最初に対応する要素 のみに範囲の設定がなされます。あるいは、変数の中の一部にのみ設定を行う場 合、下限側に -infを、上限側に infをそれぞれ lbと ub変数に設定します。た とえば、つぎの例題、 c x( )≤0 x10 ,x2≥0

(30)

lb = [-inf 0]; ub = [10 inf]; は、 に範囲の設定がなされます ( には下限がなく、 には上限 がありません)。inf および-infを使用すると、非常に大きな正の数または負の 数を使用して、範囲に制限がないことを示すよりも、より有効な数値結果を出し ます。 探索空間をさらに制限したので、解の探索のための関数の計算回数が少なくなる ことに注目してください。問題の制約および境界に関する制限が多いときには、 通常、関数値計算も少なくなります。なぜならば、最適化により、ステップサイ ズおよび可能領域に関して、制約なしの場合よりも、より良い決定を出力するか らです。そのため、解への収束を促進するには、可能な限り問題に範囲や制約を 課すのが良いと言えます。

勾配指定を伴った制約付きの例題

通常、中規模問題の最小化ルーチンでは、有限差分近似により計算した数値勾配 を使用します。この手順は、関数および制約関数の偏微分を計算するために、各 変数にシステマティックな摂動を与えます。一方、関数を与えて解析的に偏微分 を計算することができます。通常は、このような関数を設定する方が、より正確 で、かつ効率的に求めることができます。 解析的に決定した勾配を使用して、(2-2)式 を解くには、つぎのようにします。

ステップ 1: 目的関数および勾配に対する M- ファイルの記述

function [f,G] = objfungrad(x) f = exp(x(1))*(4*x(1)^2+2*x(2)^2+4*x(1)*x(2)+2*x(2)+1); % 目的関数の勾配 t = exp(x(1))*(4*x(1)^2+2*x(2)^2+4*x(1)*x(2)+2*x(2)+1); G = [ t + exp(x(1)) * (8*x(1) + 4*x(2)), exp(x(1))*(4*x(1)+4*x(2)+2)];

ステップ 2: 非線形制約および非線形制約の勾配に対する M- ファイルの

記述

function [c,ceq,DC,DCeq] = confungrad(x)

c(1) = 1.5 + x(1) * x(2) - x(1) - x(2); % 不等式制約 c(2) = -x(1) * x(2)-10; % 制約の勾配 DC= [x(2)-1, -x(2); x(1)-1, -x(1)]; x1≤10 , 0≤x2 x1 x2

(31)

標準的なアルゴリズムを使用した例題 % 非線形等式制約なし ceq=[]; DCeq = [ ]; Gは、xの各要素に対して、objfungrad(x)により求めた目的関数 fの偏微係数 を含みます。 (2-4) DCの列は、各制約に対する偏微係数を含みます ( すなわち、DCの第 i列は xに 対する第 i番目の制約の偏微係数です )。そのため、上記の例題では、DCは、つ ぎのようになります。 (2-5) objfungrad.m に目的関数の勾配、confungrad.mに制約の勾配を与えているの で、fminconに、これらの M- ファイルが追加情報を含むことを示す必要があり

ま す。optimset を 使 用 し て、パ ラ メ ー タ GradObj と GradConstr を 既 存 の

options構造体の中で、'on'に切り替えてください。

options = optimset(options,'GradObj','on','GradConstr','on');

構造体 options の中のこれらのパラメータを 'on'に設定しないと、fminconは、

解析的な勾配を使用しません。 引数 lbおよび ubは、xの中の独立変数に上下限を設定します。この例題では、境 界制約がないため、[]に設定されています。

ステップ 3: 制約付き最適化ルーチンの呼び出し

x0 = [-1,1]; % 初期値推定 options = optimset('LargeScale','off'); options = optimset(options,'GradObj','on','GradConstr','on'); fx ex1(4x12+2x22+4x1x2+2x2+1)+ex1(8x1+4x2) ex1(4x1+4x2+2) = c1x1 ∂ --- ∂c2 x1 ∂ --- c1x2 ∂ --- ∂c2 x2 ∂ --- x2–1 –x2 x1–1 –x1 =

(32)

lb = [ ]; ub = [ ]; % 上下限の設定なし [x,fval] = fmincon(@objfungrad,x0,[],[],[],[],lb,ub,... @confungrad,options) [c,ceq] = confungrad(x) % x での制約値のチェック 20 回の関数計算の後、出力される解は、つぎの通りです。 x = -9.5474 1.0474 fval = 0.0236 c= 1.0e-14 * 0.1110 -0.1776 ceq = []

勾配のチェック : 解析的に求めた勾配と数値的に求めた勾配

解析的に設定した勾配が与えられる場合、この勾配を有限差分を使って計算した ものと比較できます。これは、目的関数あるいは勾配関数の定式化のどちらに間 違いがあるのかを調べる際に特に有効です。 このような勾配チェックを行いたい場合、optimsetを使って、DerivativeCheck オプションを 'on'に設定してください。 options = optimset(options,'DerivativeCheck','on'); 最適化の最初のサイクルで、解析的に決定した ( 目的関数の ) 勾配 ( と、存在する 場合、非線形制約 ) をチェックします。設定の許容範囲内で有限差分勾配と一致 しない場合、ワーニングメッセージにより矛盾を示して、最適化の停止あるいは 続行のオプションを示します。

等式制約付き例題

等式制約の可能なルーチンでは、非線形等式制約を非線形不等式制約と共に、M-ファイル内で計算する必要があります。線形等式では、等式に関する係数として、 左辺側を行列 Aeqに、右辺側をベクトル beqに渡します。 たとえば、非線形等式制約 および非線形不等式制約 があ る場合、これらをつぎのように書き換えてください。 x12+x2 = 1 x1x2≥–10

(33)

標準的なアルゴリズムを使用した例題 そして、問題を解くには以下のようにします。

ステップ 1: M- ファイル objfun.m の記述

function f = objfun(x) f = exp(x(1))*(4*x(1)^2+2*x(2)^2+4*x(1)*x(2)+2*x(2)+1);

ステップ 2: 非線形制約に対する M- ファイル confuneq.m の記述

function [c, ceq] = confuneq(x) % 非線形不等式制約 c = -x(1)*x(2) - 10; % 非線形等式制約 ceq = x(1)^2 + x(2) - 1;

ステップ 3: 制約付き最適化ルーチンの呼び出し

x0 = [-1,1]; % 解の初期推定値 options = optimset('LargeScale','off'); [x,fval] = fmincon(@objfun,x0,[],[],[],[],[],[],... @confuneq,options) [c,ceq] = confuneq(x) % x での制約のチェック 21 回の関数計算の後、出力される解は、つぎの通りです。 x = -0.7529 0.4332 fval = 1.5093 c= -9.6739 ceq = 4.0684e-010 ceqは、制約に関するデフォルトの許容範囲 1.0e-006より小さいので 0 とみな され、cは希望したようにゼロより小さいか等しい値になっていることに注目し てください。 x12+x2–1 = 0 x1x2 – –10≤0

(34)

最大化

最 適 化 関 数 fminbnd, fminsearch, fminunc, fmincon, fgoalattain, fminimax,

lsqcurvefit, および lsqnonlinは、すべて、目的関数 の最小化を実行する も の で す。最 大 化 を 実 行 す る に は、 と し て 入 力 し ま す。 同 様 に し て、 quadprogに対する最大化を実行するには、-Hと -fを、linprogに対しては、-f を入力します。

ゼロより大きい制約

Optimization Toolbox では、非線形不等式制約の形式として、 を想定して います。0 より大きい制約を表現する場合には、両辺に -1 を乗算してください。 たとえば、形式 の制約は の制約に等しく、形式 の制約は と等価です。

ユーザ関数を入れ子関数としてパラメタライズ

ユーザ関数を無名関数としてかく代わりとして、つぎの 1 つの M- ファイルをか くことができます。 • ユーザ関数への追加のパラメータを入力としてアクセプトする • 最適化関数の呼び出し • ユーザ関数を入れ子関数として含む つぎの例題は、係数 bと c の異なる値に対して、 x3+ bx+ c の零点を見つける ために M- ファイルをかく方法を説明します。 function y = findzero(b, c, x0)

options = optimset('Display', 'off'); % 表示をオフ y = fsolve(@poly, x0, options); function y = poly(x) % 多項式の計算 y = x^3 + b*x + c; end end メイン関数 findzeroは、つぎの 2 つのことを行います。 • 多項式の零点を見つけるために関数 fzeroを呼び出します。 • fzeroによりコールされる、入れ子関数polyにおいて多項式を計算します。 f x( ) f x( ) – Ci( )x ≤0 Ci( )x ≥0 (–Ci( )x )≤0 Ci( )xb Ci( )x +b – ( )≤0

(35)

標準的なアルゴリズムを使用した例題 任意の値の係数 bと cをもつ findzeroをコールできます。その後、入れ子関数 であるため、これらは、自動的に polyに渡されます。 例として、初期の点 x0 = 0を使用して、b = 2と c = 3.5をもつ多項式の零点 を見つけるために、つぎのようにfindzeroを呼びます。 x = findzero(2, 3.5, 0) これは、零点を返します。 x = -1.0945

無名関数と入れ子関数によるグローバル変数の回避

ツールボックスの最適化関数は、つぎを含む、ユーザが定義するいくつかのタイ プの関数を使用します。 • 目的関数

• 制約関数 ( fmincon, fseminf, fgoalattain, fminimaxに対して )

• Hessian と Jacobian は、大規模fmincon, fminunc, lsqnonlin, lsqcurvefit,

fsolveに対して、それぞれ、HessMultと JacobMultを乗算します。

• 出力関数 これらの関数は、独立変数の他に追加のパラメータを必要とする場合もあります。 関数にこれらの追加のパラメータを与える方法は、2 通りあります。 • ユーザ関数をパラメタライズし、ユーザ関数をコールする無名関数に対する関 数ハンドルを作成します。これについて、 2-19ページの " ユーザ関数を無名関 数を使用してパラメタライズする " に説明します。 • ソルバをコールする外側の関数内で入れ子関数としてユーザ関数をかきます。 この方法には、 2-18ページの " ユーザ関数を入れ子関数としてパラメタライズ "に述べるように、ユーザ関数間で変数を共有できるという利点があります。

ユーザ関数を無名関数を使用してパラメタライズする

例題として、fsolve を使用して、関数 ellipj の零点を見つけたいとします。 fsolveは、目的関数が 1 つの引数をとることを期待しますが、関数ellipjは、u と mの 2 つをとります。つぎの入力により、この関数をみることができます。 type ellipj

(36)

ここでは変数 uについて解きます。mは、単に Jacobi 楕円関数を指定する第 2 の パラメータです。m = 0.5に対して、 u0 = 3の近くの零点を探すためには、ワー クスペースからmのカレント値を得る、無名関数に対する関数ハンドルを作成で きます。そのため、ソルバー fsolveがこの関数ハンドルをコールする場合、パラ メータ mを、ellipjは 2 つの引き数としてコールします。つぎのコマンドを使用 して、fsolveにこの関数ハンドルを渡します。 u0 = 3; m = 0.5;

options = optimset('Display', 'off'); % 表示をオフ x = fsolve(@(u) ellipj(u,m), u0, options)

x = 3.7081

入れ子関数を用いて変数を共有

前 の 例 題 で は、fsolve に よ り 渡 さ れ る よ り も 多 く の 引 数 を も つ 既 存 の 関 数 ellipjを使用します。ユーザ独自の関数をかいている場合、上記のテクニックを 使用できます。あるいは、入れ子関数を使用する方がより便利であることもあり ます。入れ子関数には、関数間で変数を共有できる利点があります。たとえば、目 的関数がある値を決して超えないような付加的な非線形の制約がある場合に、目 的関数を最小化したいと仮定します。入れ子関数を使用して、制約関数内で目的 関数値を再計算することを避けることができます。つぎのコードは、これを説明 します。

function [x,fval] = runsharedvalues(a,b,c,d,lower) objval = []; % 共有変数の初期化 x0 = [-1,1]; % 解の初期推定値 options = optimset('LargeScale','off'); [x, fval] = fmincon(@objfun,x0,[],[],[],[],[],[],@constrfun,options); function f = objfun(x) % 非線形目的関数

% 変数 objval は、objfun と runsharedvalues が共有します

objval = exp(x(1))*(a*x(1)^2+b*x(2)^2+c*x(1)*x(2)+d*x(2)+1); f = objval;

(37)

標準的なアルゴリズムを使用した例題

function [c,ceq] = constrfun(x) % 非線形不等式制約

% 変数 objval は、objfun と runsharedvalues が共有します c(1) = - objval + lower; c(2:3) = [1.5 + x(1)*x(2) - x(1) - x(2); -x(1)*x(2) - 10]; % 非線形等式制約 ceq = []; end end ここで、a, b, c, dを様々な値にして目的関数を実行でき、最初の非線形制約によ り、目的関数がlowerの値よりも小さくならない領域で、目的関数の最小を見つ けることが常に保証されます。 [x, fval ] = runsharedvalues(-4, 2, 4, 2, 1)

Optimization terminated: first-order optimality measure less than options.TolFun and maximum constraint violation is less than options.TolCon.

Active inequalities (to within options.TolCon = 1e-006): lower upper ineqlin ineqnonlin

1 x = -0.6507 1.3030 fval = 1 [x,fval] = runsharedvalues(4, 2, 4, 2, 0.5)

Optimization terminated: first-order optimality measure less than options.TolFun and maximum constraint violation is less than options.TolCon.

(38)

lower upper ineqlin ineqnonlin 1 x = -4.9790 1.9664 fval = 0.5000 2-32ページの "fminimax を使用する Simulink の例題 " で入れ子関数での共有変数 の他の例を見ることができます。

解析ヤコビ行列による非線形方程式

この例では、デフォルトの中規模fsolveアルゴリズムの使用方法を示します。こ こで、対象となる問題は、 • 連立非線形方程式が正方、すなわち、方程式の数と未知数の数が等しい。 • となるような 1 つの解 が存在する。 例では、fsolveを使用し、等価な連立非線形方程式を導出して解くことによりバ ナナ (Rosenbrock) 関数の最小値を求めます。 において最小値となる Rosenbrock 関数は、最適化における代表的なテスト問題です。この関数は非線形 性が高く、最急降下法を使用したとしても、きわめて緩やかにしか収束しま せん。 Rosenbrock関数は、次のように与えられます。 まず、この関数を n 次関数として一般化します。n は任意の整数の偶数値です。 この関数を、一般 Rosenbrock 関数と呼びます。一般 Rosenbrock 関数は n 個の未知 数を含む n 個の 2 乗項から成ります。 F x( ) = 0 x F x( ) = 0 f x( ) = 100 x( 2x12)2+(1–x1)2 f x( ) 100 x( 2ix2i2 1)2+(1–x2 i1)2 i=1 n 2

=

(39)

標準的なアルゴリズムを使用した例題 fsolveを使用して となるような の値を求める、つまり一般 Rosenbrock 関数の最小値を求めるには、その前に、次のような等価な連立非線形方程式とし て関数を書き換える必要があります。 この連立方程式は正方であり、fsolveを使用して解くことができます。例に示さ れるように、この連立方程式は、 で与えられる一意な解を持ち ます。

ステップ 1: 目的関数値および ヤコビ行列を計算する M- ファイル

bananaobj.m

の記述

function [F,J] = bananaobj(x); % 一般的な n 次元 Rosenbrock 関数から得られた % 非線形方程式システムのベクトル関数と % ヤコビ行列を評価する。 % 問題のサイズを取得 n = length(x);

if n == 0, error('Input vector, x, is empty.'); end if mod(n,2) ~= 0,

error('Input vector, x, must have an even number of components.'); end % ベクトル関数の計算 odds = 1:2:n; evens = 2:2:n; F = zeros(n,1); F(odds,1) = 1-x(odds); F(evens,1) = 10.*(x(evens)-x(odds).^2); F x( ) = 0 x F 1( ) = 1–x1 F 2( ) = 10 x( 2x12) F 3( ) = 1–x3 F 4( ) = 10 x( 4x32) F n( –1) = 1–xn1 F n( ) = 10 x( nxn21)

..

.

xi= 1,i=1, ,… n

(40)

% nargout >1 の場合、ヤコビ行列を計算 if nargout > 1 c = -ones(n/2,1); C = sparse(odds,odds,c,n,n); d = 10*ones(n/2,1); D = sparse(evens,evens,d,n,n); e = -20.*x(odds); E = sparse(evens,odds,e,n,n); J = C + D + E; end

ステップ 2: 方程式のシステムに対するソルバルーチンの呼び出し

n = 64; x0(1:n,1) = -1.9; x0(2:2:n,1) = 2; options=optimset('Display','iter','Jacobian','on'); [x,F,exitflag,output,JAC] = fsolve(@bananaobj,x0,options); 初期の点として、奇数インデックスに対して を使用し、偶数インデッ クスに対して を使用します。LargeScaleオプションには、fsolveのデ フォルト 'off' を使用し、中規模非線形方程式アルゴリズムにはデフォルトの 'dogleg'を使用します。そして、bananaobj.mに定義されたヤコビ行列を使用す

るために、Jacobianを 'on' に設定します。関数 fsolveは、つぎの出力をします。

Norm of First-order Trust-region Iteration Func-count f(x) step optimality radius 0 1 4281.92 615 1 1 2 1546.86 1 329 1 2 3 112.552 2.5 34.8 2.5 3 4 106.24 6.25 34.1 6.25 4 5 106.24 6.25 34.1 6.25 5 6 51.3854 1.5625 6.39 1.56 6 7 51.3854 3.90625 6.39 3.91 7 8 43.8722 0.976562 2.19 0.977 8 9 37.0713 2.44141 6.27 2.44 9 10 37.0713 2.44141 6.27 2.44 10 11 26.2485 0.610352 1.52 0.61 11 12 20.6649 1.52588 4.63 1.53 12 13 17.2558 1.52588 6.97 1.53 13 14 8.48582 1.52588 4.69 1.53 14 15 4.08398 1.52588 3.77 1.53 15 16 1.77589 1.52588 3.56 1.53 16 17 0.692381 1.52588 3.31 1.53 x i( ) = –1.9 x i( ) = 2

(41)

標準的なアルゴリズムを使用した例題

17 18 0.109777 1.16206 1.66 1.53 18 19 0 0.0468565 0 1.53 Optimization terminated successfully:

First-order optimality is less than options.TolFun

有限差分ヤコビ行列による非線形方程式

前の例題においては、関数 bananaobjが、Fを評価し、ヤコビ行列 Jを計算しま す。それでは、ヤコビ行列を計算するコードが使用できない場合はどうするので しょうか? デフォルトでは、ユーザが (optionsの中のJacobianオプションを 'on' に 設 定 し て ) 目 的 関 数 の 中 で ヤ コ ビ 行 列 を 計 算 し な い 場 合、fsolve, lsqnonlin, lsqcurvefitは有限差分を使ってヤコビ行列を近似します。これは

Jacobianオプションのデフォルトです。optimsetでJacobianオプションを'off' に設定することによって、有限差分を選択することができます。

この例では、目的関数として前の例の bananaobjを使用していますが、fsolve

がJacobianの概算値を求めbananaobjの第二出力を無視するように、Jacobian

パラメータを 'off' に設定します。LargeScale オプションには、fsolveのデ

フォルト 'off' を使用し、中規模非線形方程式アルゴリズムにはデフォルトの 'dogleg'を使用します。 n = 64; x0(1:n,1) = -1.9; x0(2:2:n,1) = 2; options=optimset('Display','iter','Jacobian','off'); [x,F,exitflag,output,JAC] = fsolve(@bananaobj,x0,options); この例題は、つぎの出力を生成します。

Norm of First-order Trust-region Iteration Func-count f(x) step optimality radius 0 65 4281.92 615 1 1 130 1546.86 1 329 1 2 195 112.552 2.5 34.8 2.5 3 260 106.24 6.25 34.1 6.25 4 261 106.24 6.25 34.1 6.25 5 326 51.3854 1.5625 6.39 1.56 6 327 51.3854 3.90625 6.39 3.91 7 392 43.8722 0.976562 2.19 0.977 8 457 37.0713 2.44141 6.27 2.44 9 458 37.0713 2.44141 6.27 2.44 10 523 26.2485 0.610352 1.52 0.61

(42)

11 588 20.6649 1.52588 4.63 1.53 12 653 17.2558 1.52588 6.97 1.53 13 718 8.48582 1.52588 4.69 1.53 14 783 4.08398 1.52588 3.77 1.53 15 848 1.77589 1.52588 3.56 1.53 16 913 0.692381 1.52588 3.31 1.53 17 978 0.109777 1.16206 1.66 1.53 18 1043 0 0.0468565 0 1.53 Optimization terminated successfully:

First-order optimality is less than options.TolFun

この例を有限差分型にした場合も、収束には前の例の解析ヤコビ行列型の場合と 同じ回数の繰り返しが必要になります。一般に、どちらの型も繰り返しにより、同 程度の割合で収束します。ただし、有限差分型では、他にも多くの関数評価処理 が必要です。このような付加的評価処理のコストが大きくなるかどうかは個々の 問題に依存します。

複数の目的関数をもつ例題

これまでの例題は、単一の目的関数を使った問題でした。本節では、lsqnonlin, fminimax, fgoalattainを使って、複数の目的関数をもつ問題の解法を示します。 最初の 2 つの例は、Simulink モデルのパラメータを最適化する方法を示します。 この節には以下の例題が含まれます。 • 2-26ページの "lsqnonlin を使用する Simulink の例題 " • 2-32ページの "fminimax を使用する Simulink の例題 " • 2-35ページの " 信号処理の例題 "

lsqnonlin

を使用する Simulink の例題

Simulink モデル optsim.mdlの制御パラメータを最適化しましょう ( このモデル

は、Optimization Toolbox のoptimディレクトリ内にあります。このモデルをロー

ドするには、Simulink がインストール済みでなければならないことに注意してく

ださい ) 。このモデルは、図 2-4, 飽和要素を含むアクチュエータをもつプラント

に示す Simulink ブロックダイアグラムとしてモデル化された非線形プロセスプ ラントを含みます。

(43)

標準的なアルゴリズムを使用した例題 図 2-4: 飽和要素を含むアクチュエータをもつプラント このプラントは、アクチュエータに制限を持つ低減衰の 3 次のモデルです。アク チュエータは、飽和制限とスルーレート制限を持ちます。アクチュエータは、入 力値を最大 2 単位、最小 -2 単位の範囲に制限しています。そして、その変化率も 0.8 単位 / 秒以内としています。ステップ入力に対するシステムの開ループ応答を 図 2-5, 閉ループ応答 に示します。モデルをオープン ( コマンドラインで optsim と入力 するか、あるいはモデル名をクリック ) し、Simulation メニューから Start を選択することにより、この応答を確認することができます。応答は Scope ブ ロックにプロット表示されます。 Actuator Model 1 y Rate 1.5 50s +a2.s +a1.s+13 2 Plant Limit 1 u

(44)

図 2-5: 閉ループ応答 問題は、システムへの単位ステップ入力をトラックするフィードバック制御ルー プの設計です。閉ループプラントは、プラントとアクチュエータをまとめて一つ の階層構造の Subsystem ブロックとして表わします。 Scope ブロックは、設計プロ セスの間、出力軌跡を表示します。図 2-6, 閉ループモデルを参照してください。 図 2-6: 閉ループモデル

(45)

標準的なアルゴリズムを使用した例題 この問題を解く 1 つの方法は、出力と入力信号間の誤差を最小化することです。 変数は PID コントローラのパラメータです。誤差の最小化をある時間単位で行う 場合、単一の目的関数を使うことができます。しかしながら、0 から 100 までの全 時間ステップでの誤差の最小化を目的とする場合、複数の目的関数 ( 各時間ス テップに対して 1 つの目的関数 ) を作成します。 ルーチン lsqnonlinを使用して、出力のトラッキングに対して最小二乗適合を実 行します。トラッキングは、M- ファイル関数 tracklsqにより実行されます。こ の関数は、sim により、計算された出力から、入力信号 1 を差し引いた誤差信号

youtを返します。下記に示す、tracklsqのコードは、 Optimization Toolboxに含

まれる、ファイルruntracklsq.mの内容です。 関数 runtracklsqは、必要な値すべてを設定し、目的関数 tracklsq を使用し て、lsqnonlinをコールします。これは、runtracklsqの中で入れ子構造になっ ています。lsqnonlinに渡す変数optionsは、種々の変数やシミュレーションパ ラメータの規範を設定したり、表示特性を定義します。出力が必要な場合には、中 規模アルゴリズムを使用してください。そして、0.001のオーダーで、ステップ や目的関数に対して終了許容範囲を設定してください。 モデルoptsimで、シミュレーションを実行するには、変数 Kp, Ki, Kd, a1, および a2 (a1と a2は Plant ブロックの中の変数 ) をすべて定義する必要があります。Kp, Ki, Kd が、最適化される変数です。関数tracklsqは、変数 a1と a2が 2 つの関 数間で共有されるように、runtracklsqの中で入れ子構造になっています。変数 a1と a2は、runtracklsq内で初期化されます。 目的関数tracklsqは、シミュレーションを実行する必要があります。シミュレー ションは、ベースワークスペースあるいはカレントのワークスぺース、すなわち、 simを呼び出している関数のワークスペース ( この場合では、tracklsqのワーク スペース)のいずれかで実行することができます。この例題では、コマンドsimset を使用して、'SrcWorkspace'を 'Current'に設定することにより、カレントの ワークスペースでシミュレーションを実行するように sim に命令します。関数 simsetを使用して、simに対するソルバを選択することもできます。このシミュ レーションは、固定ステップ 5 次法を使用して、100秒まで実行します。

シミュレーションが完了すると、変数 tout, xoutおよび xoutはカレントのワー

クスペース ( すなわち、tracklsqのワークスペース ) に存在します。ブロックダ

イアグラム内で Outport ブロックを使用して、シミュレーションの最後にyout を

カレントのワークスペースに出力します。

つぎに、runtracklsqに対するコードを示します。

function [Kp,Ki,Kd] = runtracklsq

(46)

optsim % モデルのロード

pid0 = [0.63 0.0504 1.9688]; % 初期値の設定 a1 = 3; a2 = 43; % モデル内の plant 変数の初期化

options = optimset('LargeScale','off','Display','iter',... 'TolX',0.001,'TolFun',0.001);

pid = lsqnonlin(@tracklsq, pid0, [], [], options); Kp = pid(1); Ki = pid(2); Kd = pid(3);

function F = tracklsq(pid) % 1 の信号に対する optsim の出力の追跡します % 変数 a1 と a2 は、モデル optsim により必要です % これらは、RUNTRACKLSQ と共有されているので、ここで % 再定義される必要はありません Kp = pid(1); Ki = pid(2); Kd = pid(3); % 関数値を計算 simopt = simset('solver','ode5','SrcWorkspace','Current'); % sim オプションの初期化

[tout,xout,yout] = sim('optsim',[0 100],simopt); F = yout-1; end end runtracklsq を実行する場合、最適化により、64 回の関数値計算の後、コント ローラの比例、積分、微分 (Kp, Ki, Kd)ゲインに対する解を求めます。 [Kp, Ki, Kd] = runtracklsq Directional

Iteration Func-count Residual Step-size derivative Lambda 0 4 8.66531 1 18 5.21604 85.4 -0.00836 6.92469 2 25 4.53699 1 -0.107 0.0403059 3 32 4.47316 0.973 -0.00209 0.0134348 4 40 4.46854 2.45 9.72e-005 0.00676229 5 47 4.46575 0.415 -0.00266 0.00338115 6 48 4.46526 1 -0.000999 0.00184785 Optimization terminated: directional derivative along

(47)

標準的なアルゴリズムを使用した例題

gradient less than 10*(TolFun+TolX). Kp = 3.0956 Ki = 0.1466 Kd = 14.1378 ステップ入力に対する閉ループ応答の結果を 図 2-7 に示します。 図 2-7: lsqnonlin を使用する閉ループ応答

(48)

注意 simを呼び出すことは、結果的には Simulink の常微分方程式 (ODE) ソルバ の 1 つを呼び出すこととなります。使用するソルバタイプを選択する必要があり ます。最適化という点から見ると、ODE が解けるのならば、固定ステップソルバ の選択がベストです。しかし、スティッフなシステムの場合には、可変ステップ法 により ODE を解くことが必要になることがあります。しかし、可変ステップソル バで求めた数値解は、ステップサイズ制御メカニズムのためにパラメータのス ムーズな関数ではありません。このようなスムーズ性の欠如により、最適化ルー チンが収束しなくなる場合があります。固定ステップソルバを使用すると、ス ムーズ性がなくなることはありません。( 詳細は、[1] を参照。) Simulink の可変ス テップソルバと接続して、複数の目的関数の最適化問題を解く場合、 Nonlinear

Control Design (NCD) Blockset を推奨します。これを使用すると、Simulink と共に 稼動し、特殊な数値勾配計算を行い、スムーズ性の欠如の問題を回避します。

fminimax

を使用する Simulink の例題

2-27ページの " 飽和要素を含むアクチュエータをもつプラント " に示す Simulink モデルの制御パラメータを最適化する他の方法は、関数 fminimaxを使用するこ とです。この場合、出力と入力信号間の誤差を最小にするのではなく、 0 と 100 の 間の任意の時刻tでの出力の最大値を最小にします。 下記のこの例題に対するコードは、関数 runtrackmm に含まれています。この中 で、目的関数は、単にコマンド sim により返される出力 yout です。しかしなが ら、全時間ステップで最大出力を最小化すると、幾つかの時間ステップに対して 出力が 1 よりはるかに小さくなる可能性があります。初めの20秒の後、0.95を 超える出力を保つために、制約関数 trackmmconは、t=20から t=100まで、制約 yout >= 0.95を含みます。制約は形式 g <= 0となる必要があるため、関数内で の制約は、g = -yout(20:100)+.95となります。

trackmmobjと trackmmconは、共にカレントの PID 値から計算されたsimから

の結果yout を使用します。非線形制約関数は、常に目的関数の直後に fmincon,

fminimax, fgoalattainおよび fseminfにおいて同一値とともに呼び出されま

す。こうして、入れ子関数を使用してシミュレーションを 2 度呼ぶことを回避で

き、youtの値がruntrackmmで初期化されている限り、目的関数と制約関数の間

で共有可能となります。

つぎに、runtrackmmに対するコードを示します。

function [Kp, Ki, Kd] = runtrackmm

% Copyright 1990-2004 The MathWorks, Inc.

(49)

標準的なアルゴリズムを使用した例題

optsim

pid0 = [0.63 0.0504 1.9688];

% a1, a2, yout は、TRACKMMOBJ と TRACKMMCON で共有されます a1 = 3; a2 = 43; % モデル内の plant 変数の初期化 yout = []; % yout に初期値を与える options = optimset('Display','iter',... 'TolX',0.001,'TolFun',0.001); pid = fminimax(@trackmmobj,pid0,[],[],[],[],[],[],... @trackmmcon,options);

Kp = pid(1); Ki = pid(2); Kd = pid(3); function F = trackmmobj(pid)

% 1 の信号に対する optsim の出力を追跡します % 変数 a1 と a2 は、RUNTRACKMM と共有されます

% 変数 yout は、RUNTRACKMM と RUNTRACKMMCON で共有されます Kp = pid(1);

Ki = pid(2); Kd = pid(3); % 関数値を計算

opt = simset('solver','ode5','SrcWorkspace','Current'); [tout,xout,yout] = sim('optsim',[0 100],opt);

F = yout; end

function [c,ceq] = trackmmcon(pid)

% 1 の信号に対する optsim の出力を追跡します

% 変数 yout は、 RUNTRACKMM と TRACKMMOBJ で共有されます % 制約を計算 % 目的関数 TRACKMMOBJ がこの前に呼び出されることがわかります % 制約関数、従って yout がカレントです c = -yout(20:100)+.95; ceq=[]; end end コードを実行する場合、つぎの結果が返ります。

(50)

[Kp,Ki,Kd] = runtrackmm

Max Directional

Iter F-count {F,constraints} Step-size derivative Procedure 0 5 1.11982

1 11 1.264 1 1.18 2 17 1.055 1 -0.172

3 23 1.004 1 -0.0128 Hessian modified twice 4 29 0.9997 1 3.48e-005 Hessian modified 5 35 0.9996 1 -1.36e-006 Hessian modified twice Optimization terminated: Search direction less than 2*options.TolX

and maximum constraint violation is less than options.TolCon. Active inequalities (to within options.TolCon = 1e-006): lower upper ineqlin ineqnonlin

1 14 182 Kp = 0.5894 Ki = 0.0605 Kd = 5.5295 出力の MAX{F,constraints}列に示された最後の値は、全時間ステップに対する 最大値が 0.9996 であることを表わしています。閉ループ応答の結果を 図 2-8, fminimaxを使用した閉ループ応答に示します。 この解は、問題の定式化を異なる方法で行ったlsqnonlinによる解とは異なりま す。

Figure

図 2-5:  閉ループ応答 問題は、システムへの単位ステップ入力をトラックするフィードバック制御ルー プの設計です。 閉ループプラントは、プラントとアクチュエータをまとめて一つ の階層構造の Subsystem ブロックとして表わします。  Scope  ブロックは、設計プロ セスの間、出力軌跡を表示します。 図 2-6, 閉ループモデルを参照してください。 図 2-6:  閉ループモデル

図 2-5:

閉ループ応答 問題は、システムへの単位ステップ入力をトラックするフィードバック制御ルー プの設計です。 閉ループプラントは、プラントとアクチュエータをまとめて一つ の階層構造の Subsystem ブロックとして表わします。 Scope ブロックは、設計プロ セスの間、出力軌跡を表示します。 図 2-6, 閉ループモデルを参照してください。 図 2-6: 閉ループモデル p.44
図 2-8:  fminimax を使用した閉ループ応答

図 2-8:

fminimax を使用した閉ループ応答 p.51
図 2-9:  初期および最終的な振幅係数による振幅応答00.050.10.150.20.25 0.3 0.35 0.4 0.45 0.5−3−2−10123Frequency (Hz)Magnitude Response (dB)initialfinal

図 2-9:

初期および最終的な振幅係数による振幅応答00.050.10.150.20.25 0.3 0.35 0.4 0.45 0.5−3−2−10123Frequency (Hz)Magnitude Response (dB)initialfinal p.54
図 3-3:  Rosenbrock 関数への Gauss-Newton 法の適用 Gauss-Newton  法は、 (3-16) 式の二次の項 Q(x) が意味をもつとき、しばしば問題 が生じます。 この問題を解く方法が、Levenberg-Marquardt 法です。 Levenberg-Marquardt  法 Levenberg-Marquardt  法 [27],[29] は、つぎの線形方程式の解を探索方向として使い ます。 (3-18) ここで、スカラ は、 の大きさと方向を共にコントロールし

図 3-3:

Rosenbrock 関数への Gauss-Newton 法の適用 Gauss-Newton 法は、 (3-16) 式の二次の項 Q(x) が意味をもつとき、しばしば問題 が生じます。 この問題を解く方法が、Levenberg-Marquardt 法です。 Levenberg-Marquardt 法 Levenberg-Marquardt 法 [27],[29] は、つぎの線形方程式の解を探索方向として使い ます。 (3-18) ここで、スカラ は、 の大きさと方向を共にコントロールし p.128
図 3-6:  非線形制約付き Rosenbrock 関数 (3-2) 上の SQP 法 SQP  法の実行 SQP  法は、つぎの小節で簡単に論議されている3つの部分から作られています。 • ラグランジュ関数のヘッセ行列の更新 • 二次計画問題の解法 • ライン探索とメリット関数の計算 ヘッセ行列の更新 各々の主繰り返しで、ラグランジュ行列関数のヘッセ行列の正定準 Newton 近似 H  は、 を Lagrange 乗数の推定とする BFGS 法を使って計算されま す。 -1-0.500.511.522

図 3-6:

非線形制約付き Rosenbrock 関数 (3-2) 上の SQP 法 SQP 法の実行 SQP 法は、つぎの小節で簡単に論議されている3つの部分から作られています。 • ラグランジュ関数のヘッセ行列の更新 • 二次計画問題の解法 • ライン探索とメリット関数の計算 ヘッセ行列の更新 各々の主繰り返しで、ラグランジュ行列関数のヘッセ行列の正定準 Newton 近似 H は、 を Lagrange 乗数の推定とする BFGS 法を使って計算されま す。 -1-0.500.511.522 p.138
図 3-11:   ε 制約法の幾何学的表現 このアプローチは、 重み付き和の方法を使って得ることのできない非凸境界上で、 たとえば、解  および のような、非劣位解の数を認識すること ができます。 この方法がもつ問題は、実行可能な領域解を保証するために  を適 切に選択しなければならないことです。 さらに、この方法の欠点は、厳密な制約 を設定し過ぎて、真の設計対象を表現するのにほとんど適していないことです。 似たような方法があります。 たとえば、 目的関数に優先権を設定する Waltz [40] の 方法

図 3-11:

ε 制約法の幾何学的表現 このアプローチは、 重み付き和の方法を使って得ることのできない非凸境界上で、 たとえば、解 および のような、非劣位解の数を認識すること ができます。 この方法がもつ問題は、実行可能な領域解を保証するために を適 切に選択しなければならないことです。 さらに、この方法の欠点は、厳密な制約 を設定し過ぎて、真の設計対象を表現するのにほとんど適していないことです。 似たような方法があります。 たとえば、 目的関数に優先権を設定する Waltz [40] の 方法 p.154
図 3-12:  ゴール到達法の幾何学的表示 ゴール  の指定は、ゴール点 P を定義します。 重みベクトルは、P から実 行可能関数空間 までの探索方向を定義します。 最適化の間、     は変化し、こ れ に よ り 実 行 可 能 領 域 の 大 き さ が 変 化 し ま す。 制 約 境 界 は、ユ ニ ー ク な 解 に収束します。 ゴール到達法のアルゴリズム実行 ゴール到達法は、非線形計画問題として考えられる利点をもっています。 問題の 特性は、 非線形計画アルゴリズムの中で利用されています。

図 3-12:

ゴール到達法の幾何学的表示 ゴール の指定は、ゴール点 P を定義します。 重みベクトルは、P から実 行可能関数空間 までの探索方向を定義します。 最適化の間、 は変化し、こ れ に よ り 実 行 可 能 領 域 の 大 き さ が 変 化 し ま す。 制 約 境 界 は、ユ ニ ー ク な 解 に収束します。 ゴール到達法のアルゴリズム実行 ゴール到達法は、非線形計画問題として考えられる利点をもっています。 問題の 特性は、 非線形計画アルゴリズムの中で利用されています。 p.156

References