Обеспечьте производительность LQR Используя терминальные веса штрафа

Скрипт Open Live Script

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

Стандартная функция стоимости MPC похожа на функцию стоимости для контроллера LQR с выходным взвешиванием, как показано в следующем уравнении:

$J (u) = \sum_{i = 1}^{\infty} y {(k + i)}^{T} Q y (k + i) + u {(k + i - 1)}^{T} R u (k + i - 1)$

LQR и функции стоимости MPC отличаются следующими способами:

Силы функции стоимости LQR $y$ и $u$ к нулю, тогда как функция стоимости MPC обеспечивает $y$ и $u$ к ненулевым заданным значениям. Можно переключить источник модели предсказания MPC, чтобы устранить это различие и достигнуть нулевых номинальных заданных значений.
Функция стоимости LQR использует бесконечный горизонт предсказания, в котором переменная, которой управляют, изменяется в каждом шаге расчета. В стандартной функции стоимости MPC длина горизонта $p$ , и переменные изменения, которыми управляют, $m$ времена, где $m$ горизонт управления.

Эти две функции стоимости эквивалентны, если функция стоимости MPC:

$J (u) = \sum_{i = 1}^{p - 1} (y {(k + i)}^{T} Q y (k + i) + u {(k + i - 1)}^{T} R u (k + i - 1)) + x (k + p)^{T} Q_{p} x (k + p)$

Здесь, $Q_{p}$ терминальный вес штрафа, примененный на итоговом шаге горизонта предсказания, и предсказание и горизонты управления равны ( $p$ = $m$ ). Необходимое $Q_{p}$ вычисленное использование матрицы Рикатти lqr и lqry команды.

Задайте модель объекта управления

Задайте разомкнутый цикл дискретного времени динамическая модель объекта управления с шагом расчета 0.1 секунды. Для этой модели сделайте все состояния измеримыми выходными параметрами объекта. Этот объект является объектом с двойным интегратором от [1].

A = [1 0;0.1 1];
B = [0.1;0.005];
C = eye(2);
D = zeros(2,1);
Ts = 0.1;
plant = ss(A,B,C,D,Ts);

Спроектируйте контроллер горизонта Бога LQR

Вычислите матрицу Riccati Qp и обратная связь состояния получает K сопоставленный с проблемой LQR с выходным весом Q и входной вес R. Для получения дополнительной информации смотрите lqry.

Q = eye(2);
R = 1;
[K,Qp] = lqry(plant,Q,R);

Спроектируйте эквивалентный контроллер MPC

Чтобы реализовать функцию стоимости MPC, сначала вычислите $L$ , разложение Холесского $Q_{p}$ , таким образом, что $L^{T} L = Q_{p}$ .

L = chol(Qp);

Затем задайте вспомогательные неизмеренные выходные переменные $y_{c} = L x$ , таким образом, что $y_{c}^{T} y_{c} = x^{T} Q_{p} x$ . Увеличьте выходной вектор объекта, таким образом, что это включает эти вспомогательные выходные параметры.

newPlant = plant;
set(newPlant,'C',[C;L],'D',[D;zeros(2,1)]);

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

newPlant = setmpcsignals(newPlant,'MO',[1 2],'UO',[3 4]);

Создайте объект контроллера с тем же шагом расчета как объект и равняйтесь предсказанию и управляйте горизонтами.

p = 3;
m = p;
mpcobj = mpc(newPlant,Ts,p,m);

-->The "Weights.ManipulatedVariables" property of "mpc" object is empty. Assuming default 0.00000.
-->The "Weights.ManipulatedVariablesRate" property of "mpc" object is empty. Assuming default 0.10000.
-->The "Weights.OutputVariables" property of "mpc" object is empty. Assuming default 1.00000.
   for output(s) y1 and zero weight for output(s) y2 y3 y4

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

ywt = sqrt(diag(Q))';
uwt = sqrt(diag(R))';
mpcobj.Weights.OV = [sqrt(diag(Q))' 0 0];
mpcobj.Weights.MV = sqrt(R);

Чтобы сделать проблему QP сопоставленной с контроллером MPC положительный определенный, включайте очень маленькие веса на переменном шаге, которым управляют.

mpcobj.Weights.MVRate = 1e-5;

Наложите терминальный штраф $x^{T} (k + p) Q_{p} x (k + p)$ путем определения модульного веса на $y_{c} (k + p) = L x (k + p)$ . Терминальный вес на u(t+p-1) остается то же самое.

Y = struct('Weight',[0 0 1 1]);
U = struct('Weight',uwt);
setterminal(mpcobj,Y,U);

Поскольку измеренный выходной вектор содержит целый вектор состояния, демонтируйте любой дополнительный выходной интегратор воздействия, вставленный контроллером MPC.

setoutdist(mpcobj,'model',ss(zeros(4,1)));

Удалите средство оценки состояния путем определения следующего уравнения обновления измерения:

x[n|n] = x[n|n-1] + I * (x[n]-x[n|n-1]) = x[n]

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

setEstimator(mpcobj,[],eye(2));

Сравните MPC и усиления контроллера LQR

Вычислите усиление контроллера MPC, когда ограничения будут неактивны (неограниченный MPC) и сравнят его с усилением LQR.

mpcgain = dcgain(ss(mpcobj));

-->The "Model.Noise" property of the "mpc" object is empty. Assuming white noise on each measured output channel.

fprintf('\n(unconstrained) MPC: u(k)=[%8.8g,%8.8g]*x(k)',mpcgain(1),mpcgain(2));

(unconstrained) MPC: u(k)=[-1.6355962,-0.91707456]*x(k)

fprintf('\n                LQR: u(k)=[%8.8g,%8.8g]*x(k)\n\n',-K(1),-K(2));

                LQR: u(k)=[-1.6355962,-0.91707456]*x(k)

Усиления обратной связи состояния являются точно тем же самым.

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

Сравните производительность контроллера LQR, контроллера MPC с терминальными весами и типичного контроллера MPC.

Вычислите ответ с обратной связью для контроллера LQR.

clsys = feedback(plant,K);
Tstop = 6;
x0 = [0.2;0.2];
[yLQR,tLQR] = initial(clsys,x0,Tstop);

Вычислите ответ с обратной связью для контроллера MPC с терминальными весами.

simOpt = mpcsimopt(mpcobj);
simOpt.PlantInitialState = x0;
r = zeros(1,4);
[y,t,u] = sim(mpcobj,ceil(Tstop/Ts),r,simOpt);

Создайте типичный контроллер MPC с предсказанием по умолчанию и управляйте горизонтами (p=10, m=3). Чтобы совпадать с другими контроллерами, удалите выходное возмущение и средство оценки состояния по умолчанию от типичного диспетчера MPC.

mpcobjSTD = mpc(plant,Ts);

-->The "PredictionHorizon" property of "mpc" object is empty. Trying PredictionHorizon = 10.
-->The "ControlHorizon" property of the "mpc" object is empty. Assuming 2.
-->The "Weights.ManipulatedVariables" property of "mpc" object is empty. Assuming default 0.00000.
-->The "Weights.ManipulatedVariablesRate" property of "mpc" object is empty. Assuming default 0.10000.
-->The "Weights.OutputVariables" property of "mpc" object is empty. Assuming default 1.00000.
   for output(s) y1 and zero weight for output(s) y2

mpcobjSTD.Weights.MV = uwt;
mpcobjSTD.Weights.OV = ywt;
setoutdist(mpcobjSTD,'model',tf(zeros(2,1)))
setEstimator(mpcobjSTD,[],C)

Вычислите ответ с обратной связью для типичного контроллера MPC.

simOpt = mpcsimopt(mpcobjSTD);
simOpt.PlantInitialState = x0;
r = zeros(1,2);
[ySTD,tSTD,uSTD] = sim(mpcobjSTD,ceil(Tstop/Ts),r,simOpt);

-->The "Model.Noise" property of the "mpc" object is empty. Assuming white noise on each measured output channel.

Сравните ответы контроллера.

plot(tSTD,ySTD,'r',t,y(:,1:2),'b',tLQR,yLQR,'mo')
xlabel('Time')
ylabel('Plant Outputs')
legend('Standard MPC','MPC with Terminal Weights','LQR','Location','NorthEast')

У контроллера MPC с терминальными весами есть более быстрое время урегулирования по сравнению с типичным контроллером MPC. Контроллер LQR и контроллер MPC с терминальными весами выполняют тождественно.

Можно улучшать типичный контроллер MPC производительность путем корректировки горизонтов. Например, если вы увеличиваете предсказание и управляете горизонтами (p=20, m=5), типичный контроллер MPC выполняет почти тождественно контроллеру MPC с терминальными весами.

Этот пример показывает, что использование терминальных весов штрафа может избавить от необходимости настраивать предсказание и горизонты управления для неограниченного случая MPC. Если ваше приложение включает ограничения, использование терминального веса недостаточно, чтобы гарантировать номинальную устойчивость. Необходимо также выбрать соответствующие горизонты и возможно добавить терминальные ограничения. Для получения дополнительной информации см. [2].

Ссылки

[1] Scokaert, P. O. M. и Дж. Б. Ролингс, "Ограниченное линейное квадратичное регулирование", Транзакции IEEE на Автоматическом управлении (1998), Издание 43, № 8, стр 1163-1169.

[2] Ролингс, J. B. и Д. К. Майн, прогнозирующее управление модели: теория и проект. Nob Hill Publishing, 2010.

Смотрите также

setterminal

Документация