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

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

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

この例では、強化学習 (RL) エージェントを使用して、PI コントローラーのゲインを計算します。PI コントローラーをニューラル ネットワーク コントローラーに置き換える例については、Create Simulink Environment and Train Agentを参照してください。

環境モデル

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

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)

環境インターフェイス オブジェクトを作成します。これを行うには、この例の終わりで定義されている localCreatePIDEnv を使用します。

[env,obsInfo,actInfo] = localCreatePIDEnv(mdl);

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

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

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

rng(0)

TD3 エージェントの作成

アクターを作成するには、まず観測入力とアクション出力をもつ深層ニューラル ネットワークを作成します。詳細については、rlContinuousDeterministicActorを参照してください。

PI コントローラーは、誤差および誤差積分の観測値を使用する 1 つの全結合層をもつニューラル ネットワークとしてモデル化できます。

$$\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 はアクター ニューラル ネットワークの出力。

• Kp と Ki はニューラル ネットワークの重みの絶対値。

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

勾配降下の最適化により、重みが負の値になる可能性があります。負の重みを回避するには、標準の fullyConnectedLayer を fullyConnectedPILayer に置き換えます。この層は、関数 $Y=\mathrm{abs}\left(\mathrm{WEIGHTS}\right)*X$ を実装して重みが必ず正になるようにします。これは、ファイル fullyConnectedPILayer.m (このエージェントを使用するときに MATLAB® 検索パスまたは現在のフォルダー内に存在しなければならない) で定義されています。カスタム層の定義の詳細については、カスタム深層学習層の定義を参照してください。

initialGain = single([1e-3 2]);

actorNet = [
    featureInputLayer(numObs)
    fullyConnectedPILayer(initialGain,'ActOutLyr')
    ];

actorNet = dlnetwork(actorNet);

actor = rlContinuousDeterministicActor(actorNet,obsInfo,actInfo);

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

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

クリティックを作成するには、この例の終わりで定義されている localCreateCriticNetwork を使用します。両方のクリティック表現に同じネットワーク構造を使用します。

criticNet = localCreateCriticNetwork(numObs,numAct);

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

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

critic2 = rlQValueFunction(initialize(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 を使用するようにエージェントを設定。

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

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

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

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

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

agentOpts.TargetPolicySmoothModel.StandardDeviation = sqrt(0.1);

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

agent = rlTD3Agent(actor,critic,agentOpts);

エージェントの学習

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

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

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

• 連続する 100 個を超えるエピソードで -355 よりも大きい平均累積報酬をエージェントが受け取ったときに学習を停止。この時点で、エージェントはタンクの水位を制御できます。

詳細については、rlTrainingOptionsを参照してください。

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

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

doTraining = false;

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

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

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

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

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

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

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

Ki = abs(parameters{1}(1))

Ki = single
    0.3958

Kp = abs(parameters{1}(2))

Kp = single
    8.0822

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.98024        1.7073        0.40451     10.077

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

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.9226        84.241       true  

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

rlCumulativeCost  = sum(rlCost.Data)

rlCumulativeCost = -375.9135

cstCumulativeCost = sum(cstCost.Data)

cstCumulativeCost = -376.9373

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

ローカル関数

貯水タンクの RL 環境を作成する関数。

function [env,obsInfo,actInfo] = localCreatePIDEnv(mdl)

% 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 interface object.
env = rlSimulinkEnv(mdl,[mdl '/RL Agent'],obsInfo,actInfo);

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

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

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

クリティック ネットワークを作成する関数。

function criticNet = localCreateCriticNetwork(numObs,numAct)
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');
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