Обеспечивайте эффективность LQR, используя веса штрафов терминалов

Открыть 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- Контроллера бесконечного горизонта

Вычислите матрицу Риккати 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')

Figure contains an axes. The axes contains 6 objects of type line. These objects represent Standard MPC, MPC with Terminal Weights, LQR.

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

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

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

Ссылки

[1] Scokaert, P. O. M. and J. B. Rawlings, «Ограниченное линейное квадратичное регулирование», Транзакции IEEE по автоматическому управлению (1998), том 43, № 8, стр. 1163-1169.

[2] Rawlings, J. B. and D. Q. Mayne, Model Predictive Control: Theory and Design. Nob Hill Publishing, 2010.

См. также

setterminal

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