MATLABで使える進化計算ライブラリPlatEMOの紹介（とエネルギー法による倒立振子の振り上げ）

この記事は[MATLAB/SIMULINK Advent Calendar 2024](https://qiita.com/advent-calendar/2024/matlab)の13日目の記事です！
13日の金曜日:thinking::thinking::thinking:

https://qiita.com/advent-calendar/2024/matlab

間違ってるところがあれば教えてください。。。
# 最適化
最適化とは、ある制約の下で目的関数を最小（最大）化することです。
1目的の最大化問題はこんなやつ
<img width="300" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3819721/812c1cbf-42fb-a8c8-ea23-1b7eb263a471.png">
2目的の最小化問題はこんなやつですね。
<img width="750" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3819721/4e329e9a-12b1-e064-f41c-a14393d3b7a3.png">



## 進化計算
### 遺伝的アルゴリズムとは
2目的の最小化問題を考えます。
最初は、変数$\boldsymbol{x}$を適当な値で何パターンか計算してみます。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3819721/0dac9b4f-d12d-9034-d67a-8c1d9a9fd5e2.png)

この初期個体たちは目的関数空間に散らばっていますが、最小化問題であれば左下にいる個体の方が偉く、右上にいる個体は成績の悪いダメなやつらということになります。

そして次の世代では、成績の良い個体が生き残り子孫を残し（特徴を受け継ぐ）、ダメなやつらは淘汰されます。
これを繰り返すことによって、だんだんと成績が良くなっていき、最適解へと近づいていきます。

:::note warn
局所解に陥らないように突然変異させたりというのもありますが、詳しくはググってください()
:::

# PlatEMOとは

## 概要
PlatEMOは安徽大学のBIMK (Institute of Bioinspired Intelligence and Mining Knowledge) グループによって開発された、MATLABベースの進化的多目的最適化プラットフォームです。200以上の進化アルゴリズムと400以上のベンチマーク問題を提供する、オープンソースの研究・教育用プラットフォームとなっています。

### 主な特徴
- 豊富なアルゴリズム：遺伝的アルゴリズム、差分進化、粒子群最適化など
- 充実したベンチマーク問題
- 使いやすいGUIインターフェース
- 並列実験機能
- Excel/テーブル形式での結果出力

### 技術要件
- MATLABのみで動作（追加のライブラリ不要）
- クロスプラットフォーム対応
- MATLABが動作する環境であれば利用可能

## 利用上の注意点

### ライセンスと引用
研究目的での利用は無料ですが、論文などでPlatEMOを使用した場合は以下の論文を引用する必要があります：

> Ye Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB Platform for Evolutionary Multi-Objective Optimization [Educational Forum], IEEE Computational Intelligence Magazine, 2017, 12(4): 73-87

### 最新バージョン（v4.8）の特徴
- ブロック接続による新アルゴリズムの視覚的作成
- 複数ソリューションの同時評価機能
- 新しい多目的進化アルゴリズム（MOBCA, NRV-MOEA, NSBiDiCoなど）の追加

:::note info
PlatEMOは継続的に更新されており、最新のアルゴリズムや機能が定期的に追加されています。GitHubのリポジトリをチェックすることで、最新の更新情報を確認できます。
:::

MATLABのFile Exchangeにも上がっているので、アドオンエクスプローラーからも追加できます。
![Screenshot 2024_12_12 23_03_51.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3819721/fd28220d-f87f-f6a1-da95-003ed71b6823.png)


## GUIを見てみる
PlatEMOはCUIでも実行できますが、ここではGUIでいろいろ見てみます。
![PlatEMO v4.7 2024_12_12 22_40_41.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3819721/a660c758-b3c0-13e6-a63b-7dc6d061fcea.png)
この画面では、自身で最適化問題を定義し、実行できます。
## ベンチマーク問題を遊んでみる
PlatEMOにはさまざまなベンチマーク問題があります。試しにやってみます。
ベンチマーク問題は、右上のTest Moduleからいろいろ選択できます。
![PlatEMO v4.7 2024_12_12 22_40_51.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3819721/517d545d-b6c9-857f-e0ef-5ea05b1ba742.png)
左からアルゴリズムとベンチマーク問題を選択し、目的関数と設計変数の数を決めて実行してみます。
### 2目的の最小化
![PlatEMO-v4.7-2024-12-12-22-57-18.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3819721/0204c5f6-5d65-2a16-d81f-18b1c06df1bc.gif)


### 3目的の最小化
![PlatEMO-v4.7-2024-12-12-23-00-27.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3819721/1c51a064-fc54-077b-b060-57ad159b4fe3.gif)
たまに突然変異により離れたところにプロットされる個体がいて、そのせいでパレートフロントが伸び縮みしているように見えちゃいますね。

# 自分で最適化問題を定義して試してみる
## エネルギー法による倒立振子の振り上げ
ということで、ここでは
- 制御則の切り替え角度
- LQR（最適レギュレータ）による安定化制御の重み
- エネルギー法による振り上げのゲイン
これらのパラメータを最適化によって求められないか考えてみます。

https://qiita.com/Carter/items/db4bb11ead93ea3a1d42

こちらの川田先生の記事を参考に台車型（というより、レール型？）の倒立振子の振り上げとそのアニメーションをやってみます。安定化制御と振り上げ制御、無制御の部分の切り替えだったりがうまくいかないと、いい感じに立ってくれません。

<blockquote class="twitter-tweet"><p lang="ja" dir="ltr">エネルギー法で振り子の振り上げしようとしてるけど、立ちそうで立たない<br><br>どこいくねんって感じ() <a href="https://t.co/3JKONXJMNG">pic.twitter.com/3JKONXJMNG</a></p>&mdash; ゆうきち@MATLAB Amb (@Yuki_MATLAB_Amb) <a href="https://twitter.com/Yuki_MATLAB_Amb/status/1808897976267608362?ref_src=twsrc%5Etfw">July 4, 2024</a></blockquote> <script async src="https://platform.twitter.com/widgets.js" charset="utf-8"></script>

----
### 設計変数の定義
設計変数として、以下のパラメータを最適化対象とします。


1. **LQRの重み**  
   - $q_1$: $x$（台車の位置）に対応する重み  
   - $q_2$: $\theta$（振り子の角度）に対応する重み  
   - $q_3$: $\dot{x}$（台車の速度）に対応する重み  
   - $q_4$: $\dot{\theta}$（振り子の角速度）に対応する重み  

2. **振り上げ制御のゲイン**  
   - $k_e$: エネルギー法で振り上げを行う際のゲイン  

3. **制御則の切り替え角度**  
   - $\theta_{stable}$: 安定化制御と無制御の切り替え角度  
   - $\theta_{swing}$: 無制御と振り上げ制御の切り替え角度  

これら7つのパラメータを設計変数として扱い、それぞれの範囲を以下のように定義します：

| パラメータ       | 下限 | 上限 | 
|------------------|------|------|
| $q_1$       | 1    | $10^5$  |
| $q_2$       | 1    | $10^5$  |
| $q_3$     | 1    | $10^5$  |
| $q_4$       | 1    | $10^5$  |
| $k_e$         | 0  | $10^5$   |
| $\theta_{stable}$ | 0    | 180°  |
| $\theta_{swing}$| 0    | 180°  |

---

### 目的関数の定義
次に、目的関数を定義します。最適化では、多目的最小化問題として以下を設定します：

1. **目的関数1**: 台車の位置 $x$ の2乗和  
   台車の移動距離が小さいほど評価が高くなります。
   $$
   f_1 = \sum_{t}^{t_{end}} x(t)^2
   $$

2. **目的関数2**: 振り子の角度 $\theta$ の2乗和  
   振り子が垂直に安定するほど評価が高くなります。
   $$
   f_2 = \sum_{t}^{t_{end}} \theta(t)^2
   $$



<details><summary>コードたち</summary>
ODE45でシミュレーションを行うために、事前に運動方程式などなどは求めておき、loadで呼び出しています。

:::note warn
ここでは2つの目的関数両方でODE45を実行しているので無駄が多いです。一回で済ませるにはグローバル変数とかキャッシュとか使うんですかね？
:::

``` matlab:costFunction.m
function [cost] = costFunction(variables)
    q1   = variables(1);
    q2   = variables(2);
    q3   = variables(3);
    q4   = variables(4);
    k_e  = variables(5);
    theta_stable = variables(6);
    theta_swing = variables(7);

    load('workspaceData.mat','A','B','StateEqu','params','x0');
    
    % 状態重み行列と入力重み
    Q = diag([q1, q2, q3, q4]);
    R = 1; 

    % LQRゲイン計算
    K = lqr(A, B, Q, R);

    % シミュレーション設定
    tf = 10;
    tspan = 0:0.01:tf;

    % ode45を用いて非線形システムシミュレーション
    [t, x] = ode45(@(t,x) pendulum_cart_nonlinear(t, x, K, StateEqu, params, k_e, theta_stable, theta_swing), tspan, x0);

    % コスト計算 xの2乗和
    cost = sumsqr(x(:,1));

end
```


``` matlab:costFunction2.m
function [cost] = costFunction2(variables)
    q1   = variables(1);
    q2   = variables(2);
    q3   = variables(3);
    q4   = variables(4);
    k_e  = variables(5);
    theta_stable = variables(6);
    theta_swing = variables(7);

    load('workspaceData.mat','A','B','StateEqu','params','x0');
    
    % 状態重み行列と入力重み
    Q = diag([q1, q2, q3, q4]);
    R = 1; 

    % LQRゲイン計算
    K = lqr(A, B, Q, R);

    % シミュレーション設定
    tf = 10;
    tspan = 0:0.01:tf;

    % ode45を用いて非線形システムシミュレーション
    [t, x] = ode45(@(t,x) pendulum_cart_nonlinear(t, x, K, StateEqu, params, k_e, theta_stable, theta_swing), tspan, x0);

    % コスト計算 thetaの2乗和
    cost = sumsqr(x(:,2));

end
```

</details>



---
### 制約条件
制約条件として、入力する力の大きさに上限を設定します。


### 初期個体
自分でもいろいろとパラメータを調整していたところ、なんとか振り上げできるパラメータを発見したので初期個体の一つにこれを与えてみます。

![InvPenAnim.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3819721/4bb329c6-f24f-54db-0a5b-131a9c233c2f.gif)

``` matlab:myInitialization.m
function PopDec = myInitialization(N)

    goodSolution = [1,10,1,1000,10,30,135];

    D = 7;
    PopDec = zeros(N,D);

    % 最初の1個体に良好解を割り当て
    PopDec(1,:) = goodSolution;

    % 残りのN-1個はランダムに生成する
    Lower = [1,1,1,1,0,0,0];
    Upper = [1e5,1e5,1e5,1e5,1e5,180,180];
    PopDec(2:end,:) = rand(N-1,D).*(Upper - Lower) + Lower;
end
```

---
### シミュレーションと最適化の実行
NSGA-IIという多目的最適化用のアルゴリズムを選択し、最適化を実行してみます。
### 結果
100個体*100世代の10000回の評価で、最終世代ではこのような結果になりました。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3819721/df9c9540-90c9-7d92-be2d-adb5d0c6dffb.png)

f2の値は全体でみるとほとんど変わっておらず、f1の値の方が個体ごとに大きく違うように見えます。f1の小さい個体の動画を見てみると、いい感じに最適化されてるなーという感じです。
![f2_opti.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3819721/ff105fc8-69de-9e81-3140-a255fce2d513.gif)

# まとめ
急いで書いたのでコードの詳細など掘り切れていません。。。CUIで使用する場合、使い方がplatemo.mの中にも書いてあるのでそちらを参照してみてください！！！

<details><summary>platemo.m</summary>

``` matlab platemo.m
function varargout = platemo(varargin)
%platemo - The main function of PlatEMO.
%
%   platemo() displays the GUI of PlatEMO.
%
%   platemo('Name',Value,'Name',Value,...) runs an algorithm on a problem
%   with specified parameter settings.
%
% All the acceptable names and values are:
%	'algorithm'     <function handle>	an algorithm
%	'problem'       <function handle>	a problem
%   'N'             <positive integer>  population size
%   'M'             <positive integer>  number of objectives
%   'D'             <positive integer>  number of variables
%	'maxFE'         <positive integer>  maximum number of function evaluations
%   'maxRuntime'    <positive integer>  maximum runtime (in second)
%   'save'       	<integer>           number of saved populations
%   'run'           <integer>           current run number
%   'metName'       <string>            names of metrics to calculate
%   'outputFcn'     <function handle>   function called after each iteration
%   'encoding'      <string>            encoding scheme of each decision variable (1.real 2.integer 3.label 4.binary 5.permutation)
%   'lower'         <vector>            lower bound of each decision variable
%   'upper'         <vector>            upper bound of each decision variable
%   'initFcn'       <function handle>   function for initializing solutions
%   'evalFcn'       <function handle>   function for evaluating solutions
%   'decFcn'        <function handle>   function for repairing invalid solutions
%   'objFcn'        <function handle>   objective functions
%   'conFcn'        <function handle>   constraint functions
%   'objGradFcn'    <function handle>   function for calculating the gradients of objectives
%   'conGradFcn'    <function handle>   function for calculating the gradients of constraints
%
%   Example:
%
%       platemo()
%
%   displays the GUI of PlatEMO.
%
%       platemo('algorithm',@GA,'problem',@SOP_F1,'N',50,'maxFE',20000)
%
%   runs GA with a population size of 50 on SOP_F1 for 20000 evaluations.
%
%       platemo('algorithm',@PSO,'problem',@SOP_F1,'N',100,'maxRuntime',3)
%
%   runs PSO with a population size of 100 on SOP_F1 for 3 seconds.
%
%       platemo('algorithm',{@KnEA,0.4},'problem',{@WFG4,6},'M',5)
%
%   runs KnEA on 5-objective WFG4 and sets the parameters in KnEA and WFG4.
%
%       for i = 1 : 10
%           platemo('algorithm',@MOEAD,'problem',@ZDT1,'save',5,'run',i)
%       
```
</details>