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

Открыть сценарий в реальном времени

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

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

$J (u)_{}^{} {=∑i=1∞y (k}^{} + i) TQy {(k + i)}^{} + u (k +$ i-1) TRu (k + i-1)

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

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

Две функции затрат эквивалентны, если функция затрат MPC:

$J (u)_{}^{} =∑i=1p-1 {(y (}^{k} + i) TQy {(k + i}^{)} + u (k + i-1)^{TRu}_{(} k +$ i-1)) + x (k + p) TQpx (k + p)

Здесь $_{Qp}$ является конечным штрафным весом, применяемым на этапе конечного горизонта прогнозирования, и горизонты прогнозирования и управления равны ( $p$ = $m$ ). Требуемый $_{Qp}$ представляет собой матрицу Риккати, рассчитанную с использованием 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$ , разложение по Чолескому $_{Qp}$ , так что $^{} LTL =_{}$ Qp.

L = chol(Qp);

Затем определите вспомогательные неизмеренные выходные переменные $_{yc} =$ Lx, так ${что}_{}^{}_{} {ycTyc}^{} =_{}$ xTQpx. Увеличьте выходной вектор установки таким образом, чтобы он включал эти вспомогательные выходы.

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;

Наложить штраф на терминал $^{xT} (k +_{p}) Qpx$ (k + p) путем указания $_{} веса$ единицы измерения для yc (k + p) = Lx (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] Скокаерт, П. О. М. и Дж. Б. Ролингс, «Ограниченная линейная квадратичная регуляция», IEEE Transactions on Automatic Control (1998), том 43, № 8, стр. 1163-1169.

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

См. также

setterminal

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