強化学習を使用した PI コントローラーの調整

この例では、双生遅延深層決定論的方策勾配 (TD3) 強化学習アルゴリズムを使用して、PI コントローラーの 2 つのゲインを調整する方法を示します。調整したコントローラーのパフォーマンスは、"制御システム調整器" アプリを使用して調整したコントローラーのパフォーマンスと比較されます。"制御システム調整器" アプリを使用して Simulink® でコントローラーを調整するには、Simulink Control Design™ ソフトウェアが必要です。

調整可能なパラメーターの数が少ない比較的単純な制御タスクの場合、モデルベースの調整手法は、モデルフリーの RL ベースの手法と比較して、より高速な調整プロセスで良好な結果を得ることができます。一方、RL 手法は、非線形性の高いシステムや適応コントローラーの調整により適している可能性があります。コントローラーの比較を容易にするために、どちらの調整方法でも線形 2 次ガウシアン (LQG) 目的関数を使用します。DDPG エージェントを使用して LQR コントローラーを実装する例については、Compare DDPG Agent to LQR Controllerを参照してください。

この例では、強化学習 (RL) エージェントを使用して、PI コントローラーのゲインを計算します。PI コントローラーをニューラル ネットワーク コントローラーに置き換える例については、DDPG エージェントを使用したタンク内の水位の制御を参照してください。

再現性のための乱数ストリームの固定

サンプル コードには、さまざまな段階で乱数の計算が含まれる場合があります。サンプル コード内のさまざまなセクションの先頭にある乱数ストリームを固定すると、実行するたびにセクション内の乱数シーケンスが維持されるため、結果が再現される可能性が高くなります。詳細については、Results Reproducibilityを参照してください。

シード 0 で乱数ストリームを固定し、メルセンヌ・ツイスター乱数アルゴリズムを使用します。乱数生成に使用されるシード制御の詳細については、rngを参照してください。

previousRngState = rng(0,"twister");

出力 previousRngState は、ストリームの前の状態に関する情報が格納された構造体です。例の最後で、状態を復元します。

環境モデル

この例の環境モデルは貯水タンク モデルです。この制御システムの目的は、基準値と一致するようにタンク内の水位を維持することです。

open_system('watertankLQG')

モデルには、分散が $\mathit{E}\left({\mathit{n}}^2\left(\mathit{t}\right)\right)=1$ のプロセス ノイズが含まれています。

制御量 u を最小限に抑えながら水位を維持するために、この例のコントローラーは次の LQG 基準を使用します。

$$\mathit{J}=\underset{\mathit{T}\Rightarrow \infty }{\lim }\mathit{E}\left(\frac{1}{\mathit{T}}{\int }_0^{\mathit{T}}\left({\left(\mathit{Href}-\mathit{y}\left(\mathit{t}\right)\right)}^2+0.01{\mathit{u}}^2\left(\mathit{t}\right)\right)\mathrm{dt}\right)$$

このモデルでコントローラーをシミュレーションするには、シミュレーション時間 Tf とコントローラーのサンプル時間 Ts を秒単位で指定しなければなりません。

Ts = 0.1;
Tf = 10;

貯水タンク モデルの詳細については、watertank Simulink モデル (Simulink Control Design)を参照してください。

制御システム調整器を使用した PI コントローラーの調整

"制御システム調整器" を使用して Simulink でコントローラーを調整するには、コントローラー ブロックを調整済みブロックとして指定し、調整プロセスの目標を定義しなければなりません。"制御システム調整器" の使用の詳細については、制御システム調整器を使用した制御システムの調整 (Simulink Control Design)を参照してください。

この例では、"制御システム調整器" を使用して、保存済みのセッション ControlSystemTunerSession.mat を開きます。このセッションでは、watertankLQG モデルの PID Controller ブロックを調整済みのブロックとして指定しており、LQG 調整目標を含んでいます。

controlSystemTuner("ControlSystemTunerSession")

コントローラーを調整するには、[調整] タブで [調整] をクリックします。

調整された比例ゲインと積分ゲインは、それぞれ約 9.8 と 1e-6 です。

Kp_CST = 9.80199999804512;
Ki_CST = 1.00019996230706e-06;

エージェント学習用の環境の作成

RL エージェントに学習させるためのモデルを定義するには、次の手順を使用して貯水タンク モデルを変更します。

1. PID Controller を削除します。

2. RL Agent ブロックを挿入します。

3. 観測値ベクトル ${\left[\begin{array}{ccc}\int \mathit{e}\mathrm{dt}& \mathit{e}& \end{array}\right]}^{\mathit{T}}$ を作成します。ここで、$\mathit{e}=\mathit{Href}-\mathit{h}$ であり、$\mathit{h}$ はタンク内の水位、$\mathit{Href}$ は基準水位です。観測信号を RL Agent ブロックに接続します。

4. RL エージェントの報酬関数を LQG コストの "負"、つまり $\mathrm{Reward}=-\left({\left(\mathit{Href}-\mathit{h}\left(\mathit{t}\right)\right)}^2+0.01{\mathit{u}}^2\left(\mathit{t}\right)\right)$ として定義します。RL エージェントはこの報酬を最大化することによって LQG コストを最小化します。

結果のモデルは rlwatertankPIDTune.slx です。

mdl = 'rlwatertankPIDTune';
open_system(mdl)

rlSimulinkEnv を使用して、貯水タンクの RL 環境オブジェクトを作成します。環境のリセット関数を設定するには、この例の終わりで定義されている localResetFcn を使用します。

% Define the observation specification obsInfo 
% and the action specification actInfo.
obsInfo = rlNumericSpec([2 1]);
obsInfo.Name = 'observations';
obsInfo.Description = 'integrated error and error';

actInfo = rlNumericSpec([1 1]);
actInfo.Name = 'PID output';

% Build the environment object.
env = rlSimulinkEnv(mdl,[mdl '/RL Agent'],obsInfo,actInfo);

% Set a custom reset function that randomizes 
% the reference values for the model.
env.ResetFcn = @(in)localResetFcn(in,mdl);

この環境の観測値およびアクションの次元を抽出します。prod(obsInfo.Dimension) と prod(actInfo.Dimension) を使用して、それぞれ観測空間と行動空間の次元の数 (行ベクトル、列ベクトル、行列のいずれによって構成されているかにかかわらず) を取得します。

numObs = prod(obsInfo.Dimension);
numAct = prod(actInfo.Dimension);

再現性をもたせるために、乱数発生器のシードを固定します。

rng(0)

TD3 エージェントの作成

再現性のために乱数ストリームを固定します。

rng(0,"twister");

PI コントローラーをアクターとしてモデル化できます。アクターは、誤差積分と誤差観測値を入力として使用します。

$$\mathit{u}=\left[\begin{array}{ccc}\int \mathit{e}\mathrm{dt}& \mathit{e}& \end{array}\right]*{\left[\begin{array}{ccc}{\mathit{K}}_{\mathit{i}}& {\mathit{K}}_{\mathit{p}}& \end{array}\right]}^{\mathit{T}}$$

ここで、以下となります。

• u はアクターの出力。

• Ki と Kp はアクターの重み。

• $\mathit{e}=\mathit{Href}-\mathit{h}\left(\mathit{t}\right)$。$\mathit{h}\left(\mathit{t}\right)$ はタンク内の水位、$\mathit{Href}$ は基準水位。

アクターを作成するには、例のフォルダーに別個のファイルとして用意されている基底関数 myBasisFcn と共に、rlContinuousDeterministicActorを使用します。詳細については、rlContinuousDeterministicActorを参照してください。

initialGain = single([1e-3 2]);
actor = rlContinuousDeterministicActor( ...
    {@myBasisFcn,initialGain'},obsInfo, actInfo);

この例のエージェントは、双生遅延深層決定論的方策勾配 (TD3) エージェントです。TD3 エージェントは、アクター近似器オブジェクトとクリティック近似器オブジェクトによって最適な方策を学習します。

TD3 エージェントは、2 つのクリティック価値関数表現を使用し、観測値とアクションに基づいて長期報酬を近似します。クリティックを作成するには、まず観測値とアクションの 2 つの入力および 1 つの出力をもつ深層ニューラル ネットワークを作成します。

クリティック ネットワークを作成します。両方のクリティック表現に同じネットワーク構造を使用します。

statePath = [
    featureInputLayer(numObs,Name='stateInLyr')
    fullyConnectedLayer(32,Name='fc1')
    ];

actionPath = [
    featureInputLayer(numAct,Name='actionInLyr')
    fullyConnectedLayer(32,Name='fc2')
    ];

commonPath = [
    concatenationLayer(1,2,Name='concat')
    reluLayer
    fullyConnectedLayer(32)
    reluLayer
    fullyConnectedLayer(1,Name='qvalOutLyr')
    ];

criticNet = dlnetwork();
criticNet = addLayers(criticNet,statePath);
criticNet = addLayers(criticNet,actionPath);
criticNet = addLayers(criticNet,commonPath);

criticNet = connectLayers(criticNet,'fc1','concat/in1');
criticNet = connectLayers(criticNet,'fc2','concat/in2');

指定したニューラル ネットワーク、および環境のアクション仕様と観測仕様を使用して、クリティック オブジェクトを作成します。追加の引数として、観測チャネルとアクション チャネルに接続するネットワーク層の名前も渡します。

critic1 = rlQValueFunction(criticNet, ...
    obsInfo,actInfo,...
    ObservationInputNames='stateInLyr', ...
    ActionInputNames='actionInLyr');

critic2 = rlQValueFunction(criticNet, ...
    obsInfo,actInfo,...
    ObservationInputNames='stateInLyr', ...
    ActionInputNames='actionInLyr');

クリティック オブジェクトから成るベクトルを定義します。

critic = [critic1 critic2];

アクターとクリティックの学習オプションを指定します。

actorOpts = rlOptimizerOptions( ...
    LearnRate=1e-3, ...
    GradientThreshold=1);

criticOpts = rlOptimizerOptions( ...
    LearnRate=1e-3, ...
    GradientThreshold=1);

rlTD3AgentOptionsを使用して TD3 エージェントのオプションを指定します。アクターとクリティックの学習オプションを含めます。

• コントローラーのサンプル時間 Ts を使用するようにエージェントを設定。

• ミニバッチ サイズを 256 経験サンプルに設定。

• 経験バッファーの長さを 1e6 に設定。

• 分散 0.1 のガウス ノイズを使用するように探索モデルとターゲット方策平滑化モデルを設定。

agentOpts = rlTD3AgentOptions(...
    SampleTime=Ts,...
    MiniBatchSize=256, ...
    ExperienceBufferLength=1e6,...
    ActorOptimizerOptions=actorOpts,...
    CriticOptimizerOptions=criticOpts);

ドット表記を使用してエージェントのオプションを設定または変更することもできます。

agentOpts.TargetPolicySmoothModel.StandardDeviation = sqrt(0.1);

指定したアクター表現、クリティック表現、およびエージェントのオプションを使用して TD3 エージェントを作成します。詳細については、rlTD3AgentOptionsを参照してください。

agent = rlTD3Agent(actor,critic,agentOpts);

エージェントの学習

再現性のために乱数ストリームを固定します。

rng(0,"twister");

エージェントに学習させるには、まず次の学習オプションを指定します。

• 最大 1000 個のエピソードについて、各学習を実行 (各エピソードの持続時間は最大 100 タイム ステップ)。

• [強化学習の学習モニター] で学習の進行状況を表示し (Plots オプションを設定)、コマンド ラインの表示を無効化 (Verbose オプションを設定)。

• 50 回の学習エピソードごとに貪欲方策のパフォーマンスを評価し、5 回のシミュレーションの累積報酬を平均化。

• 評価スコアが -345 に達したら学習を停止。この時点で、エージェントはタンクの水位を制御できます。

学習オプションの詳細については、rlTrainingOptionsを参照してください。

maxepisodes = 1000;
maxsteps = ceil(Tf/Ts);
trainOpts = rlTrainingOptions(...
    MaxEpisodes=maxepisodes, ...
    MaxStepsPerEpisode=maxsteps, ...
    ScoreAveragingWindowLength=100, ...
    Verbose=false, ...
    Plots="training-progress",...
    StopTrainingCriteria="EvaluationStatistic",...
    StopTrainingValue=-345);

% agent evaluator
evaluator = rlEvaluator( ...
    NumEpisodes=5,EvaluationFrequency=50,RandomSeeds=1:5);

関数 train を使用して、エージェントに学習させます。このエージェントの学習は計算量が多いプロセスのため、完了するのに数分かかります。この例の実行時間を節約するために、doTraining を false に設定して事前学習済みのエージェントを読み込みます。エージェントに学習させるには、doTraining を true に設定します。

doTraining = false;
if doTraining
    % Train the agent.
    trainingStats = train(agent,env,trainOpts,Evaluator=evaluator);
else
    % Load pretrained agent for the example.
    load("WaterTankPIDtd3.mat","agent")
end

学習済みエージェントの検証

再現性のために乱数ストリームを固定します。

rng(0,"twister");

シミュレーションを実行し、学習済みエージェントをモデルに対して検証します。

simOpts = rlSimulationOptions(MaxSteps=maxsteps);
experiences = sim(env,agent,simOpts);

PI コントローラーの積分ゲインと比例ゲインは、アクター表現の重みです。重みを取得するには、まずアクターから学習可能なパラメーターを抽出します。

actor = getActor(agent);
parameters = getLearnableParameters(actor);

コントローラーのゲインを取得します。

Ki = parameters{1}(1)

Ki = single

0.4862

Kp = parameters{1}(2)

Kp = single

8.2851

RL エージェントから取得したゲインを元の PI コントローラー ブロックに適用し、ステップ応答シミュレーションを実行します。

mdlTest = 'watertankLQG';
open_system(mdlTest); 
set_param([mdlTest '/PID Controller'],'P',num2str(Kp))
set_param([mdlTest '/PID Controller'],'I',num2str(Ki))
sim(mdlTest)

シミュレーションのステップ応答情報、LQG コスト、安定余裕を抽出します。安定余裕を計算するには、この例の終わりで定義されている localStabilityAnalysis を使用します。

rlStep = simout;
rlCost = cost;
rlStabilityMargin = localStabilityAnalysis(mdlTest);

"制御システム調整器" を使用して取得したゲインを元の PI コントローラー ブロックに適用し、ステップ応答シミュレーションを実行します。

set_param([mdlTest '/PID Controller'],'P',num2str(Kp_CST))
set_param([mdlTest '/PID Controller'],'I',num2str(Ki_CST))
sim(mdlTest)
cstStep = simout;
cstCost = cost;
cstStabilityMargin = localStabilityAnalysis(mdlTest);

コントローラーのパフォーマンスの比較

各システムのステップ応答をプロットします。

figure
plot(cstStep)
hold on
plot(rlStep)
grid on
legend('Control System Tuner','RL',Location="southeast")
title('Step Response')

両方のシミュレーションのステップ応答を解析します。

rlStepInfo = stepinfo(rlStep.Data,rlStep.Time);
cstStepInfo = stepinfo(cstStep.Data,cstStep.Time);
stepInfoTable = struct2table([cstStepInfo rlStepInfo]);
stepInfoTable = removevars(stepInfoTable,{'SettlingMin', ...
    'TransientTime','SettlingMax','Undershoot','PeakTime'});
stepInfoTable.Properties.RowNames = {'CST','RL'};
stepInfoTable

stepInfoTable=2×4 table
           RiseTime    SettlingTime    Overshoot     Peak 
           ________    ____________    _________    ______

    CST    0.77737        1.3278        0.33125     9.9023
    RL     0.95173        1.6188        0.45174     10.103

両方のシミュレーションの安定性を解析します。

stabilityMarginTable = struct2table( ...
    [cstStabilityMargin rlStabilityMargin]);
stabilityMarginTable = removevars(stabilityMarginTable,{...
    'GMFrequency','PMFrequency','DelayMargin','DMFrequency'});
stabilityMarginTable.Properties.RowNames = {'CST','RL'};
stabilityMarginTable

stabilityMarginTable=2×3 table
           GainMargin    PhaseMargin    Stable
           __________    ___________    ______

    CST      8.1616        84.123       true  
    RL       9.6843        83.828       true  

2 つのコントローラーの累積 LQG コストを比較します。RL 調整済みのコントローラーが、より最適な解を生成しています。

rlCumulativeCost  = sum(rlCost.Data)

rlCumulativeCost = 
-375.4192

cstCumulativeCost = sum(cstCost.Data)

cstCumulativeCost = 
-376.9373

どちらのコントローラーも安定した応答を生成しており、"制御システム調整器" を使用して調整したコントローラーはより高速な応答を生成します。一方、RL 調整手法ではより高いゲイン余裕が得られ、より最適な解が得られます。

previousRngState に保存されている情報を使用して、乱数ストリームを復元します。

rng(previousRngState);

ローカル関数

各エピソードの開始時に基準信号とタンク内の初期水位をランダム化する関数。

function in = localResetFcn(in,mdl)

% Randomize reference signal
blk = sprintf([mdl '/Desired \nWater Level']);
hRef = 10 + 4*(rand-0.5);
in = setBlockParameter(in,blk,'Value',num2str(hRef));

% Randomize initial water height
hInit = rand;
blk = [mdl '/Water-Tank System/H'];
in = setBlockParameter(in,blk,'InitialCondition',num2str(hInit));

end

SISO 貯水タンク システムの安定余裕を線形化して計算する関数。

function margin = localStabilityAnalysis(mdl)

io(1) = linio([mdl '/Sum1'],1,'input');
io(2) = linio([mdl '/Water-Tank System'],1,'openoutput');
op = operpoint(mdl);
op.Time = 5;
linsys = linearize(mdl,io,op);

margin = allmargin(linsys);
end

参考

関数

• train | sim | rlSimulinkEnv

オブジェクト

• rlTD3Agent | rlTD3AgentOptions | rlTrainingOptions

トピック

• Compare DDPG Agent to LQR Controller
• Train TD3 Agent for PMSM Control
• Create and Train Custom LQR Agent
• 双生遅延深層決定論的 (TD3) 方策勾配エージェント
• カスタム深層学習層の定義
• Train Reinforcement Learning Agents