Сгенерируйте премиальную функцию от прогнозирующего диспетчера модели для сервопривода

В этом примере показано, как автоматически сгенерировать премиальную функцию от стоимости и ограничительных технических требований, заданных в прогнозирующем объекте контроллера модели. Вы затем используете сгенерированную премиальную функцию, чтобы обучить агента обучения с подкреплением.

Введение

Можно использовать generateRewardFunction чтобы сгенерировать вознаграждение функционируют для обучения с подкреплением, начинающего со стоимости и ограничений, заданных в прогнозирующем контроллере модели. Получившимся сигналом вознаграждения является сумма затрат (как задано целевой функцией) и штрафы нарушения ограничений в зависимости от текущего состояния среды.

Этот пример основан на theDC Сервоприводе с Ограничением на Неизмеренный Выход (Model Predictive Control Toolbox) пример, в котором вы проектируете прогнозирующий контроллер модели для сервомеханизма DC под напряжением и ограничениями крутящего момента вала. Здесь, вы преобразуете стоимость и ограничительные технические требования, заданные в mpc объект в премиальную функцию и использование это, чтобы обучить агента управлять сервоприводом.

Откройте модель Simulink для этого примера, который является на основе вышеупомянутого примера MPC, но был изменен для обучения с подкреплением.

open_system('rl_motor')

Создайте прогнозирующий контроллер модели

Создайте динамическую модель разомкнутого контура двигателя, заданного в plant и максимальный допустимый крутящий момент tau использование функции помощника.

[plant,tau] = mpcmotormodel;

Задайте типы сигнала ввода и вывода для контроллера MPC. Угловое положение вала, измеряется, как сначала выведено. Второй выход, крутящий момент, неизмерим.

plant = setmpcsignals(plant,'MV',1,'MO',1,'UO',2);

Задайте ограничения на переменную, которой управляют и задайте масштабный коэффициент.

MV = struct('Min',-220,'Max',220,'ScaleFactor',440);

Наложите ограничения крутящего момента во время первых трех шагов предсказания и задайте масштабный коэффициент и для положения вала и для крутящего момента.

OV = struct('Min',{-Inf, [-tau;-tau;-tau;-Inf]},...
            'Max',{Inf, [tau;tau;tau;Inf]},...
            'ScaleFactor',{2*pi, 2*tau});

Задайте веса для квадратичной функции стоимости, чтобы достигнуть отслеживания углового положения. Обнулите вес для крутящего момента, таким образом, позволив ему плавать в рамках его ограничения.

Weights = struct('MV',0,'MVRate',0.1,'OV',[0.1 0]);

Создайте контроллер MPC для plant модель с шагом расчета 0.1 s, горизонт предсказания 10 шаги и горизонт управления 2 шаги, с помощью ранее заданных структур для весов, управляли переменными и выходными переменными.

mpcobj = mpc(plant,0.1,10,2,Weights,MV,OV);

Отобразите технические требования контроллера.

mpcobj
 
MPC object (created on 25-Aug-2021 20:00:40):
---------------------------------------------
Sampling time:      0.1 (seconds)
Prediction Horizon: 10
Control Horizon:    2

Plant Model:        
                                      --------------
      1  manipulated variable(s)   -->|  4 states  |
                                      |            |-->  1 measured output(s)
      0  measured disturbance(s)   -->|  1 inputs  |
                                      |            |-->  1 unmeasured output(s)
      0  unmeasured disturbance(s) -->|  2 outputs |
                                      --------------
Indices:
  (input vector)    Manipulated variables: [1 ]
  (output vector)        Measured outputs: [1 ]
                       Unmeasured outputs: [2 ]

Disturbance and Noise Models:
        Output disturbance model: default (type "getoutdist(mpcobj)" for details)
         Measurement noise model: default (unity gain after scaling)

Weights:
        ManipulatedVariables: 0
    ManipulatedVariablesRate: 0.1000
             OutputVariables: [0.1000 0]
                         ECR: 10000

State Estimation:  Default Kalman Filter (type "getEstimator(mpcobj)" for details)

Constraints:
 -220 <= MV1 (V) <= 220, MV1/rate (V) is unconstrained,       MO1 (rad) is unconstrained
                                                        -78.54 <= UO1 (Nm)(t+1) <= 78.54
                                                        -78.54 <= UO1 (Nm)(t+2) <= 78.54
                                                        -78.54 <= UO1 (Nm)(t+3) <= 78.54
                                                          UO1 (Nm)(t+4) is unconstrained

Контроллер работает с объектом с 4 состояниями, 1 вход (напряжение) и 2 выходных сигнала (угол и крутящий момент) и имеет следующие технические требования:

  • Веса функции стоимости для переменного, плавающего курса, которым управляют, которым управляют, и выходных переменных 0, 0.1 и [0.1 0] соответственно.

  • Переменная, которой управляют, ограничивается между-220V и 220 В.

  • Плавающий курс, которым управляют, неограничен.

  • Первая выходная переменная (угол) является неограниченной, но второй (крутящий момент), ограничивается между-78.54 нм и 78,54 нм в первых трех временных шагах предсказания и неограниченный на четвертом шаге.

Обратите внимание на то, что для обучения с подкреплением только ограничительная спецификация от первого временного шага предсказания будет использоваться, поскольку вознаграждение вычисляется для одного временного шага.

Сгенерируйте премиальную функцию

Сгенерируйте премиальный код функции из технических требований в mpc объект с помощью generateRewardFunction. Код отображен в редакторе MATLAB.

generateRewardFunction(mpcobj)

Сгенерированная премиальная функция является начальной точкой для премиального проекта. Можно изменить функцию с различным выбором функции штрафа и настроить веса. В данном примере внесите следующее изменение в сгенерированный код:

  • Масштабируйте исходные веса стоимости Qy, Qmv и Qmvrate на коэффициент 100.

  • Внешним методом функции штрафа по умолчанию является step. Измените метод в quadratic.

После того, как вы внесете изменения, технические требования стоимости и штрафа должны быть следующие:

Qy      = [10 0];
Qmv     = 0;
Qmvrate = 10;
Py      = Wy      * exteriorPenalty(y,ymin,ymax,'quadratic');
Pmv     = Wmv     * exteriorPenalty(mv,mvmin,mvmax,'quadratic');
Pmvrate = Wmvrate * exteriorPenalty(mv-lastmv,mvratemin,mvratemax,'quadratic');

В данном примере модифицированный код был сохранен в файле функции MATLAB rewardFunctionMpc.m. Отобразите сгенерированную премиальную функцию.

type rewardFunctionMpc.m
function reward = rewardFunctionMpc(y,refy,mv,refmv,lastmv)
% REWARDFUNCTIONMPC generates rewards from MPC specifications.
%
% Description of input arguments:
%
% y : Output variable from plant at step k+1
% refy : Reference output variable at step k+1
% mv : Manipulated variable at step k
% refmv : Reference manipulated variable at step k
% lastmv : Manipulated variable at step k-1
%
% Limitations (MPC and NLMPC):
%     - Reward computed based on first step in prediction horizon.
%       Therefore, signal previewing and control horizon settings are ignored.
%     - Online cost and constraint update is not supported.
%     - Custom cost and constraint specifications are not considered.
%     - Time varying cost weights and constraints are not supported.
%     - Mixed constraint specifications are not considered (for the MPC case).

% Reinforcement Learning Toolbox
% 02-Jun-2021 16:05:41

%#codegen

%% Specifications from MPC object
% Standard linear bounds as specified in 'States', 'OutputVariables', and
% 'ManipulatedVariables' properties
ymin = [-Inf -78.5398163397448];
ymax = [Inf 78.5398163397448];
mvmin = -220;
mvmax = 220;
mvratemin = -Inf;
mvratemax = Inf;

% Scale factors as specified in 'States', 'OutputVariables', and
% 'ManipulatedVariables' properties
Sy  = [6.28318530717959 157.07963267949];
Smv = 440;

% Standard cost weights as specified in 'Weights' property
Qy      = [0.1 0];
Qmv     = 0;
Qmvrate = 0.1;

%% Compute cost
dy      = (refy(:)-y(:)) ./ Sy';
dmv     = (refmv(:)-mv(:)) ./ Smv';
dmvrate = (mv(:)-lastmv(:)) ./ Smv';
Jy      = dy'      * diag(Qy.^2)      * dy;
Jmv     = dmv'     * diag(Qmv.^2)     * dmv;
Jmvrate = dmvrate' * diag(Qmvrate.^2) * dmvrate;
Cost    = Jy + Jmv + Jmvrate;

%% Penalty function weight (specify nonnegative)
Wy = [1 1];
Wmv = 10;
Wmvrate = 10;

%% Compute penalty
% Penalty is computed for violation of linear bound constraints.
%
% To compute exterior bound penalty, use the exteriorPenalty function and
% specify the penalty method as 'step' or 'quadratic'.
%
% Alternatively, use the hyperbolicPenalty or barrierPenalty function for
% computing hyperbolic and barrier penalties.
%
% For more information, see help for these functions.
%
% Set Pmv value to 0 if the RL agent action specification has
% appropriate 'LowerLimit' and 'UpperLimit' values.
Py      = Wy      * exteriorPenalty(y,ymin,ymax,'step');
Pmv     = Wmv     * exteriorPenalty(mv,mvmin,mvmax,'step');
Pmvrate = Wmvrate * exteriorPenalty(mv-lastmv,mvratemin,mvratemax,'step');
Penalty = Py + Pmv + Pmvrate;

%% Compute reward
reward = -(Cost + Penalty);
end

Чтобы интегрировать эту премиальную функцию, откройте блок MATLAB function в модели Simulink.

open_system('rl_motor/Reward Function')

Добавьте функцию со следующей строкой кода и сохраните модель.

r = rewardFunctionMpc(y,refy,mv,refmv,lastmv);

Блок MATLAB function теперь выполнит rewardFunctionMpc.m в процессе моделирования.

В данном примере блок MATLAB function был уже изменен и сохранен.

Создайте среду обучения с подкреплением

Динамика окружения моделируется в подсистеме Сервомеханизма. Для этой среды,

  • Наблюдения являются ссылочными и фактическими выходными переменными (угол и крутящий момент) от последних 8 временных шагов.

  • Действие является напряжением V примененный сервопривод.

  • Шаг расчета Ts=0.1s.

  • Общее время симуляции Tf=20s.

Задайте общее время симуляции и шаг расчета.

Tf = 20;
Ts = 0.1;

Создайте спецификации наблюдений и спецификации действия для среды.

numObs = 32;
numAct = 1;
oinfo = rlNumericSpec([numObs 1]);
ainfo = rlNumericSpec([numAct 1],'LowerLimit',-220,'UpperLimit',220);

Создайте среду обучения с подкреплением с помощью rlSimulinkEnv функция.

env = rlSimulinkEnv('rl_motor','rl_motor/RL Agent',oinfo,ainfo);

Создайте агента обучения с подкреплением

Зафиксируйте случайный seed для воспроизводимости.

rng(0)

Агент в этом примере является Двойным Задержанным Глубоким Детерминированным Градиентом политики (TD3) агент.

Создайте два представления критика.

% Critic
cnet = [
    featureInputLayer(numObs,'Normalization','none','Name', 'State')
    fullyConnectedLayer(128, 'Name', 'fc1')
    concatenationLayer(1,2,'Name','concat')
    reluLayer('Name','relu1')
    fullyConnectedLayer(64, 'Name', 'fc3')
    reluLayer('Name','relu2')
    fullyConnectedLayer(1, 'Name', 'CriticOutput')]; 
actionPath = [
    featureInputLayer(numAct,'Normalization','none', 'Name', 'Action')
    fullyConnectedLayer(8, 'Name', 'fc2')];
criticNetwork = layerGraph(cnet);
criticNetwork = addLayers(criticNetwork, actionPath);
criticNetwork = connectLayers(criticNetwork,'fc2','concat/in2');
criticOptions = rlRepresentationOptions('LearnRate',1e-3,'GradientThreshold',1);
critic1 = rlQValueRepresentation(criticNetwork,oinfo,ainfo,...
    'Observation',{'State'},'Action',{'Action'},criticOptions);
critic2 = rlQValueRepresentation(criticNetwork,oinfo,ainfo,...
    'Observation',{'State'},'Action',{'Action'},criticOptions);

Создайте представление актера.

actorNetwork = [
    featureInputLayer(numObs,'Normalization','none','Name','State')
    fullyConnectedLayer(128, 'Name','actorFC1')
    reluLayer('Name','relu1')
    fullyConnectedLayer(64, 'Name','actorFC2')
    reluLayer('Name','relu2')
    fullyConnectedLayer(numAct,'Name','Action')];
actorOptions = rlRepresentationOptions('LearnRate',1e-3,'GradientThreshold',1);
actor = rlDeterministicActorRepresentation(actorNetwork,oinfo,ainfo,...
    'Observation',{'State'},'Action',{'Action'},actorOptions);

Задайте опции агента с помощью rlTD3AgentOptions. Агент обучает от буфера опыта максимальной способности 1e6 путем случайного выбора мини-пакетов размера 256. Коэффициент дисконтирования 0.995 способствует долгосрочным вознаграждениям.

agentOpts = rlTD3AgentOptions("SampleTime",Ts, ...
    "DiscountFactor", 0.995, ...
    "ExperienceBufferLength",1e6, ...
    "MiniBatchSize",256);

Модель исследования в этом агенте TD3 является Гауссовой. Шумовая модель добавляет универсальное случайное значение в действие во время обучения. Установите стандартное отклонение шума к 100. Стандартное отклонение затухает по курсу 1e-5 каждый шаг агента до минимального значения 0.005.

agentOpts.ExplorationModel.StandardDeviationMin = 0.005;
agentOpts.ExplorationModel.StandardDeviation = 100;
agentOpts.ExplorationModel.StandardDeviationDecayRate = 1e-5;

Создайте агента TD3 с помощью представлений актёра и критика. Для получения дополнительной информации об агентах TD3 смотрите rlTD3Agent.

agent = rlTD3Agent(actor,[critic1,critic2],agentOpts);

Обучите агента

Чтобы обучить агента, сначала задайте опции обучения с помощью rlTrainingOptions. В данном примере используйте следующие опции:

  • Запустите каждое обучение самое большее 2 000 эпизодов с каждым эпизодом, длящимся в большей части ceil(Tf/Ts) временные шаги.

  • Остановите обучение, когда агент получит среднее совокупное вознаграждение, больше, чем -2 по 20 последовательные эпизоды. На данном этапе агент может отследить опорный сигнал.

trainOpts = rlTrainingOptions(...
    'MaxEpisodes',2000, ...
    'MaxStepsPerEpisode',ceil(Tf/Ts), ...
    'StopTrainingCriteria','AverageReward',...
    'StopTrainingValue',-2,...
    'ScoreAveragingWindowLength',20);

Обучите агента с помощью train функция. Обучение этот агент является в вычислительном отношении интенсивным процессом, который может занять несколько минут, чтобы завершиться. Чтобы сэкономить время при выполнении этого примера, загрузите предварительно обученного агента установкой doTraining к false. Чтобы обучить агента самостоятельно, установите doTraining к true.

doTraining = false;
if doTraining
    trainingStats = train(agent,env,trainOpts);
else
    load('rlDCServomotorTD3Agent.mat')
end

Снимок состояния процесса обучения показывают в следующем рисунке. Можно ожидать различные результаты из-за свойственной случайности в учебном процессе.

Подтвердите ответ контроллера

Чтобы подтвердить производительность обученного агента, симулируйте модель и просмотрите ответ в блоках Scope. Агент обучения с подкреплением может отследить опорный угол при удовлетворении ограничениям на крутящий момент и напряжение.

sim('rl_motor');

Закройте модель.

close_system('rl_motor')

Copyright 2021 The MathWorks, Inc..