Моделируйте прогнозирующее управление объектом с одним выходом с мультивходами
Этот пример показывает, как спроектировать прогнозирующий контроллер модели для объекта с одним измеренным выходом, одной манипулированной переменной, одним измеренным нарушением порядка и одним неизмеренным нарушением порядка.
Определите модель объекта управления
Определите модель объекта управления. В этом примере используются непрерывные передаточные функции времени от каждого входа к выходу. Отобразите модель объекта управления в командной строке.
plantTF = tf({1,1,1},{[1 .5 1],[1 1],[.7 .5 1]}) % display transfer functions
plantTF =
From input 1 to output:
1
---------------
s^2 + 0.5 s + 1
From input 2 to output:
1
-----
s + 1
From input 3 to output:
1
-------------------
0.7 s^2 + 0.5 s + 1
Continuous-time transfer function.
В данном примере явным образом преобразуйте объект в форму дискретного состояния пространства перед переходом к функции создания контроллера MPC.
Функция создания контроллера может принимать объекты в непрерывном времени или дискретном времени. Когда настраивается задача оптимизации, чтобы найти оптимальное значение манипулируемой переменной, контроллер MPC автоматически преобразует объект непрерывного времени (в любом формате) в модель дискретного состояния пространства, используя Нуль Order Hold.
Преобразование объекта в дискретное время полезно, когда вам нужны системные матрицы дискретного времени для анализа или симуляции (как в этом примере) или хотите, чтобы контроллер использовал объект дискретного времени, преобразованный с методом, отличным от ZOH.
plantCSS = ss(plantTF); % convert plant from transfer function to continuous-time state space Ts = 0.2; % specify a sample time of 0.2 seconds plantDSS = c2d(plantCSS,Ts); % convert plant to discrete-time state space, using Zero Order Hold
По умолчанию все входные сигналы объекта приняты как переменные с манипуляциями. Использование setmpcsignal определить, какие сигналы измеряются и не измеряются нарушения порядка. В этом примере первый входной сигнал является манипулированной переменной, второй - измеренным нарушением порядка, третий - неизмеренным нарушением порядка. Эта информация хранится в модели объекта управления plantDSS и позже используется контроллером mpc.
plantDSS = setmpcsignals(plantDSS,'MV',1,'MD',2,'UD',3); % specify signal types
Проектирование контроллера MPC
Создайте объект контроллера, задав время дискретизации, а также горизонты предсказания и управления (10 и 3 шага соответственно).
mpcobj = mpc(plantDSS,Ts,10,3);
-->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.
Поскольку вы не задали веса квадратичной функции затрат, которая будет минимизирована контроллером, их значение принято как значение по умолчанию (0 для манипулируемых переменных, 0,1 для манипулируемых переменных скоростей, 1 для выходных переменных). Кроме того, на данной точке задача MPC все еще без ограничений, так как вы еще не указали никаких ограничений.
Задайте жесткие ограничения для управляемой переменной.
mpcobj.MV = struct('Min',0,'Max',1,'RateMin',-10,'RateMax',10);
Модели входного и выходного возмущения определяют динамические характеристики неизмеренных нарушений порядка на входе и выходе, соответственно, поэтому они могут быть лучше отклонены. По умолчанию эти модели возмущения приняты как интеграторы, если вы не задаете их иначе. Объект mpc также имеет модель шума, которая задает динамику шума на измеренной выходной переменной. По умолчанию это принято как статическое усиление единицы, что эквивалентно принятию, что шум белый. В этом примере нет измеренного выходного нарушения порядка, поэтому нет необходимости задавать выходную модель возмущения, и шум на измеренном выходном сигнале принимается белым.
Задайте модель возмущения для неизмеренного входа в качестве интегратора, управляемого белым шумом с отклонением = 1000.
mpcobj.Model.Disturbance = tf(sqrt(1000),[1 0]);
Отобразите объект контроллера MPC mpcobj чтобы просмотреть его свойства.
mpcobj
MPC object (created on 27-Jan-2021 07:33:56):
---------------------------------------------
Sampling time: 0.2 (seconds)
Prediction Horizon: 10
Control Horizon: 3
Plant Model:
--------------
1 manipulated variable(s) -->| 5 states |
| |--> 1 measured output(s)
1 measured disturbance(s) -->| 3 inputs |
| |--> 0 unmeasured output(s)
1 unmeasured disturbance(s) -->| 1 outputs |
--------------
Indices:
(input vector) Manipulated variables: [1 ]
Measured disturbances: [2 ]
Unmeasured disturbances: [3 ]
(output vector) Measured outputs: [1 ]
Disturbance and Noise Models:
Output disturbance model: default (type "getoutdist(mpcobj)" for details)
Input disturbance model: user specified (type "getindist(mpcobj)" for more details)
Measurement noise model: default (unity gain after scaling)
Weights:
ManipulatedVariables: 0
ManipulatedVariablesRate: 0.1000
OutputVariables: 1
ECR: 100000
State Estimation: Default Kalman Filter (type "getEstimator(mpcobj)" for details)
Constraints:
0 <= MV1 <= 1, -10 <= MV1/rate <= 10, MO1 is unconstrained
Отобразите входную модель возмущения. Как и ожидалось, указанный интегратор преобразуется в дискретное время.
getindist(mpcobj)
ans =
A =
x1
x1 1
B =
Noise#1
x1 0.8
C =
x1
UD1 7.906
D =
Noise#1
UD1 0
Sample time: 0.2 seconds
Discrete-time state-space model.
Отобразите выходную модель возмущения. Как и ожидалось, он пуст.
getoutdist(mpcobj)
Assuming no disturbance added to measured output channel #1. -->The "Model.Noise" property of the "mpc" object is empty. Assuming white noise on each measured output channel. ans = Empty state-space model.
Исследуйте статическое смещение
Чтобы выяснить, сможет ли контроллер MPC отклонить постоянные выходные нарушения порядка и отследить постоянную уставку с нулевыми смещениями в установившемся состоянии, вычислите коэффициент усиления постоянного тока замкнутого цикла от выходных нарушений порядка до управляемых выходов с помощью cloffset команда. Это также известно как устойчивое состояние выхода чувствительность закрытого цикла.
DC = cloffset(mpcobj);
fprintf('DC gain from output disturbance to output = %5.8f (=%g) \n',DC,DC);
DC gain from output disturbance to output = 0.00000000 (=3.21965e-15)
Нулевой коэффициент усиления (который обычно является результатом того, что контроллер имеет интеграторы в качестве входных или выходных моделей возмущения) означает, что измеренный выход объекта управления будет отслеживать желаемую выходную контрольную уставку.
Симулируйте реакцию замкнутой системы с помощью sim Команда
The sim команда предоставляет быстрый способ моделирования контроллера MPC в замкнутом цикле с линейным инвариантным по времени объектом, когда ограничения и веса остаются постоянными, и нарушения порядка и опорные сигналы могут быть легко и полностью заданы раньше времени.
Во-первых, задайте время симуляции и ссылку и сигналы нарушения порядка
Tstop = 30; % simulation time Nf = round(Tstop/Ts); % number of simulation steps r = ones(Nf,1); % output reference signal v = [zeros(Nf/3,1);ones(2*Nf/3,1)]; % measured input disturbance signal
Запустите симуляцию замкнутой системы и постройте график результатов. Объект, указанный в mpcobj.Model.Plant используется как модель объекта управления в замкнутом цикле, так и как внутренняя модель объекта управления, используемая контроллером для предсказания отклика над горизонтом прогноза. Использование sim для симуляции системы замкнутого цикла для шагов Nf с ссылка и измеренного входного нарушения порядка v.
sim(mpcobj,Nf,r,v) % simulate plant and controller in closed loop
Манипулируемая переменная первоначально попадает на верхнюю границу и приводит выход объекта к ссылке значению в течение нескольких секунд. Затем он оседает на своем максимально допустимом значении, 1. Через 10 секунд измеренный нарушением порядка сигнал повышается с 0 до 1, что заставляет выход объекта превысить ссылку значения примерно на 30%. Манипулируемая переменная попадает на нижнюю границу в попытке отклонить нарушение порядка. Контроллер способен вернуть выход объекта к ссылке значению через несколько секунд, и управляемая переменная оседает на своем минимальном значении. Неизмеренный сигнал нарушения порядка всегда равен нулю, потому что неизмеренное нарушение порядка еще не задано.
Можно использовать объекты опций симуляции, чтобы задать дополнительные опции симуляции, а также сигналы шума и нарушения порядка, которые поступают в объект, но неизвестны контроллеру. В этом примере используйте объект опции симуляции, чтобы задать неизмеренный входной сигнал нарушения порядка, а также аддитивный шум как на манипулируемых переменных, так и на измеренных выходах. Вы также будете позже использовать этот объект опций, чтобы задать объект для использования в симуляции, который отличается от используемого внутри контроллера. Создайте объект опции симуляции с опциями по умолчанию и без сигнала.
SimOptions = mpcsimopt; % create object
Создайте сигнал нарушения порядка и задайте его в объекте опций симуляции
d = [zeros(2*Nf/3,1);-0.5*ones(Nf/3,1)]; % define a step disturbance signal SimOptions.UnmeasuredDisturbance = d; % specify unmeasured input disturbance signal
Задайте сигналы шума в объекте опций симуляции. Во времени симуляции функция симуляции непосредственно добавляет заданный выходной шум к измеренному выходу перед подачей его на контроллер. Он также непосредственно добавляет заданный входной шум к манипулируемой переменной (не к никаким сигналам нарушения порядка) перед подачей его на объект.
SimOptions.OutputNoise=.001*(rand(Nf,1)-.5); % specify output measurement noise SimOptions.InputNoise=.05*(rand(Nf,1)-.5); % specify noise on manipulated variables
Обратите внимание, что вы также можете использовать OutputNoise поле объекта опции симуляции для определения более общего дополнительного выходного сигнала нарушения порядка (такого как, например, шаг) на выходе измеряемой установки.
Использование sim с дополнительными SimOptions аргумент для симуляции системы замкнутого цикла и сохранения результатов в переменных рабочей области y, |t|,|u| и xp. Это позволяет вам избирательно отображать сигналы в новом окне рисунка и в любом заданном цвете и порядке.
[y,t,u,xp] = sim(mpcobj,Nf,r,v,SimOptions); % simulate closed loop
Постройте график результатов.
figure % create new figure subplot(2,1,1) % create upper subplot plot(0:Nf-1,y,0:Nf-1,r) % plot plant output and reference title('Output') % add title so upper subplot ylabel('MO1') % add a label to the upper y axis grid % add a grid to upper subplot subplot(2,1,2) % create lower subplot plot(0:Nf-1,u) % plot manipulated variable title('Input'); % add title so lower subplot xlabel('Simulation steps') % add a label to the lower x axis ylabel('MV1') % add a label to the lower y axis grid % add a grid to lower subplot
Несмотря на добавленный шум (который особенно заметен на манипулируемом переменном графике) и несмотря на измеренные и немеренные нарушения порядка, вводимые после 50 и 100 шагов соответственно, контроллер способен достигать и поддерживать хорошее отслеживание. Манипулируемая переменная оседает примерно в 1 после начальной части симуляции (стадии от 20 до 50), примерно в 0 для отклонения измеренного нарушения порядка (стадии от 70 до 100) и примерно в 0,5 для отклонения обоих нарушений порядка (стадии от 120 до 150).
Симулируйте реакцию замкнутой системы с несоответствием модели
Проверьте робастность контроллера MPC на несоответствие модели. Задайте истинное объект для использования в симуляции следующим truePlantCSS.
truePlantTF = tf({1,1,1},{[1 .8 1],[1 2],[.6 .6 1]}) % specify and display transfer functions
truePlantCSS = ss(truePlantTF); % convert to continuous state space
truePlantCSS = setmpcsignals(truePlantCSS,'MV',1,'MD',2,'UD',3); % specify signal types
truePlantTF =
From input 1 to output:
1
---------------
s^2 + 0.8 s + 1
From input 2 to output:
1
-----
s + 2
From input 3 to output:
1
-------------------
0.6 s^2 + 0.6 s + 1
Continuous-time transfer function.
Обновите объект опции симуляции путем определения SimOptions.Model как структура с двумя полями, Plant (содержащая истинную модель объекта управления) и Nominal (содержащее значения рабочих точек для истинного объекта). В данном примере номинальные значения для производных по состоянию и входов не заданы, поэтому они приняты равными нулю, что приводит к y = SimOptions.Model.Nominal.Y + C*(x-SimOptions.Model.Nominal.X) где x и y - состояние и измеренный выход объекта управления, соответственно.
SimOptions.Model = struct('Plant',truePlantCSS); % create the structure and assign the 'Plant' field SimOptions.Model.Nominal.Y = 0.1; % create and assign the 'Nominal.Y' field SimOptions.Model.Nominal.X = -.1*[1 1 1 1 1]; % create and assign the 'Nominal.X' field SimOptions.PlantInitialState = [0.1 0 -0.1 0 .05]; % specify the initial state of the true plant
удалите любой сигнал, ранее добавленный к измеренному выходу и к манипулируемой переменной
SimOptions.OutputNoise = []; % remove output measurement noise SimOptions.InputNoise = []; % remove noise on manipulated variable
Запустите симуляцию замкнутой системы и постройте график результатов. Начиная с SimOptions.Model не пуст, SimOptions.Model.Plant преобразуется в дискретное время (с использованием удержания нулевого порядка) и используется в качестве модели объекта управления в симуляции замкнутого цикла, в то время как объект в mpcobj.Model.Plant используется контроллером только для предсказания отклика над горизонтом предсказания.
sim(mpcobj,Nf,r,v,SimOptions) % simulate the closed loop
-->Converting model to discrete time.
В результате несоответствия модели видна некоторая деградация в отклике, особенно контроллеру нужно немного больше времени, чтобы достичь отслеживания, и манипулируемая переменная теперь оседает на уровне около 0,5, чтобы отклонить измеренное нарушение порядка (см. значения от 5 до 10 секунд) и устанавливается на уровне около 0,9, чтобы отклонить обоих входных нарушений порядка (от 25 до 30 секунд). Несмотря на это, контроллер по-прежнему способен отслеживать выход ссылку.
Смягчение ограничений
Каждое ограничение связано с параметром без размерности, называемым значением ECR. Ограничение с большим значением ECR допускается нарушать больше, чем ограничение с меньшим значением ECR. По умолчанию все ограничения, накладываемые на управляемые переменные, имеют значение ECR, равное нулю, что делает их жесткими. Можно задать ненулевое значение ECR для ограничения, чтобы сделать его мягким.
Обратите внимание, что можно обратиться к ManipulatedVariables также как MV. Ослабьте ограничения на управляемые переменные от жестких до мягких.
mpcobj.MV.MinECR = 1; % ECR for the lower bound on the manipulated variable mpcobj.MV.MaxECR = 1; % ECR for the upped bound on the manipulated variable
Задайте ограничение выхода. По умолчанию все ограничения, накладываемые на выходные переменные (измеренные выходы), имеют значение ECR, равное единице, что делает их мягкими. Можно уменьшить значение ECR для ограничения выхода, чтобы сделать его более трудным, однако всегда рекомендуется сохранять выходные ограничения мягкими. Это происходит главным образом из-за того, что выходы объекта зависят от состояния объекта, а также от измеренных нарушений порядка, поэтому, если происходит большое нарушение порядка, выходы объекта могут нарушать ограничения независимо от любого управляющего действия, предпринятого контроллером mpc, особенно когда управляемые переменные сами по себе жестко ограничены. Такое неизбежное нарушение жесткого ограничения приводит к недопустимой задаче MPC, для которой никакое MV не может быть вычислено.
mpcobj.OV.Max = 1.1; % define the (soft) output constraint % Note that is possible to refer to the |OutputVariables| field also as |OV|.
Запустите новую симуляцию с обратной связью, не включая объект опции симуляции, и, следовательно, без какого-либо несоответствия модели, неизмеренного нарушения порядка или шума, добавленного к манипулируемой переменной или измеренному выходу.
sim(mpcobj,Nf,r,v) % simulate the closed loop
Assuming no disturbance added to measured output channel #1. -->The "Model.Noise" property of the "mpc" object is empty. Assuming white noise on each measured output channel.
В попытке отклонить измеренное нарушение порядка, добиться отслеживания и предотвратить увеличение выходов выше своего мягкого ограничения 1.1, контроллер слегка нарушает мягкое ограничение на манипулируемую переменную, которое достигает небольших отрицательных значений с 10 до 11 секунд (вы можете увеличить изображение, чтобы увидеть это нарушение более четко). Ограничение на измеряемом выходе нарушается больше, чем ограничение на манипулируемой переменной.
Штрафуйте больше выход нарушений ограничений и перезапустите симуляцию.
mpcobj.OV.MaxECR = 0.001; % the closer to zero, the harder the constraint sim(mpcobj,Nf,r,v) % run a new closed-loop simulation.
Assuming no disturbance added to measured output channel #1. -->The "Model.Noise" property of the "mpc" object is empty. Assuming white noise on each measured output channel.
Теперь контроллер нарушает ограничение выхода только очень слабо. Как ожидалось, это улучшение эффективности на выходе происходит за счет нарушения манипулируемого ограничения переменной намного больше (манипулируемая переменная достигает -3 на пару шагов).
Изменение коэффициентов усиления Калмана, используемых во встроенной оценке состояния
На каждом временном шаге контроллер MPC получает управляемую переменную путем умножения состояния объекта в дискретном времени (при наличии) на матрицу усиления (вычисленную путем решения ограниченной квадратичной задачи оптимизации). Поскольку состояние объекта обычно недоступно, по умолчанию контроллер использует линейный фильтр Калмана в качестве наблюдателя, чтобы оценить состояние объекта с увеличенным дискретным временем (то есть объекта, комплексного с моделями нарушения порядка и шума). Поэтому состояниями контроллера являются состояния этого фильтра Калмана, которые в свою очередь являются оценками состояний дополненного объекта дискретного времени.
Запустите симуляцию замкнутого цикла с несовпадением модели и неизмеренным нарушением порядка, используя оценщик по умолчанию, и верните структуру состояний контроллера xc в рабочей области, чтобы можно было позже сравнить последовательности состояний.
[y,t,u,xp,xc] = sim(mpcobj,Nf,r,v,SimOptions);
-->Converting model to discrete time.
отобразить компонент структуры состояний контроллера
xc
xc =
struct with fields:
Plant: [150x5 double]
Disturbance: [150x1 double]
Noise: [150x0 double]
LastMove: [150x1 double]
Постройте график выхода ответа объекта управления, а также состояний объекта, оцененных наблюдателем по умолчанию
figure; % create figure subplot(2,1,1) % create upper subplot axis plot(t,y) % plot y versus time title('Plant output'); % add title to upper plot ylabel('y') % add a label to the upper y axis grid % add grid to upper plot subplot(2,1,2) % create lower subplot axis plot(t,xc.Plant) % plot xc.Plant versus time title('Estimated plant states'); % add title to lower plot xlabel('Time (sec)') % add a label to the lower x axis ylabel('xc') % add a label to the lower y axis grid % add grid to lower plot
Как ожидалось, происходят внезапные изменения на 10 и 20 секундах из-за того, что измеренные и не измеренные сигналы нарушения порядка вводятся, соответственно.
Можно изменить усиления фильтра Калмана, используемого в качестве наблюдателя. Для этого сначала найдите стандартные матрицы усиления Калмана и пространства состояний.
[L,M,A1,Cm1] = getEstimator(mpcobj); % retrieve observer matrices
Вычислите и отобразите полюса наблюдателя по умолчанию. Обратите внимание, что, как и ожидалось, все они находятся внутри модуля круга, хотя некоторые из них кажутся относительно близкими к границе. Обратите внимание, что существует шесть состояний, пять принадлежащих модели объекта управления, шестое - модели входа возмущения.
e = eig(A1-A1*M*Cm1) % eigenvalues of observer state matrix
e = 0.5708 + 0.4144i 0.5708 - 0.4144i 0.4967 + 0.0000i 0.9334 + 0.1831i 0.9334 - 0.1831i 0.8189 + 0.0000i
Спроектируйте новую оценку состояния путем размещения полюса.
poles = [.8 .75 .7 .85 .6 .81]; % specify desired positions for the new poles L = place(A1',Cm1',poles)'; % calculate Kalman gain for time update M = A1\L; % calculate Kalman gain for measurement update
Установите коэффициент усиления новой матрицы в объекте контроллера mpc
setEstimator(mpcobj,L,M); % set the new estimation gains
повторный запуск симуляции замкнутого цикла
[y,t,u,xp,xc] = sim(mpcobj,Nf,r,v,SimOptions);
Assuming no disturbance added to measured output channel #1. -->The "Model.Noise" property of the "mpc" object is empty. Assuming white noise on each measured output channel. -->Converting model to discrete time.
Постройте график выхода ответа объекта управления, а также состояний объекта, оцененных новым наблюдателем
figure; % create figure subplot(2,1,1) % create upper subplot axis plot(t,y) % plot y versus time title('Plant output'); % add title to upper plot ylabel('y') % add a label to the upper y axis grid % add grid to upper plot subplot(2,1,2) % create lower subplot axis plot(t,xc.Plant) % plot xc.Plant versus time title('Estimated plant states'); % add title to lower plot xlabel('Time (sec)') % add a label to the lower x axis ylabel('xc') % add a label to the lower y axis grid % add grid to lower plot
Как ожидалось, состояния контроллера отличаются от таковых, нанесенных ранее, и общая реакция закрытого цикла несколько медленнее.
Симулируйте ответ разомкнутого контура
Проверяйте поведение объекта и контроллера в разомкнутом контуре, используя sim команда. Установите OpenLoop флаг в on, и обеспечивают последовательность манипулируемых переменных, чтобы возбуждать систему.
SimOptions.OpenLoop = 'on'; % set open loop option SimOptions.MVSignal = sin((0:Nf-1)'/10); % define the manipulated variable signal
Симулируйте систему разомкнутого контура, содержащую истинное объект (ранее указанный в SimOptions.Model), за которым следует контроллер. Поскольку опорный сигнал будет проигнорирован, можно избежать его указания.
sim(mpcobj,Nf,[],v,SimOptions) % simulate the open loop system
-->Converting model to discrete time.
Симулируйте контроллер в замкнутом цикле с помощью mpcmove
В более общем случае закрытого цикла, в котором контроллер применяется к нелинейному или изменяющемуся во времени объекту, ограничения или веса могут измениться во время исполнения, или нарушению порядка и опорные сигналы трудно задать полностью раньше времени, можно использовать mpcmove функция на каждом шаге в цикле for для симуляции контроллера MPC в закрытом цикле. Если ваш объект находится в непрерывном режиме, вам нужно будет преобразовать его в дискретное время (с использованием любого метода), чтобы симулировать его с циклом for.
Сначала получите матрицы пространства состояний в дискретном времени объекта управления и задайте время симуляции и начальные состояния для объекта и контроллера.
[A,B,C,D] = ssdata(plantDSS); % discrete-time plant plant ss matrices Tstop = 30; % Simulation time x = [0 0 0 0 0]'; % Initial state of the plant xmpc = mpcstate(mpcobj); % Initial state of the MPC controller r = 1; % Output reference signal
Отображение начального состояния контроллера. Обратите внимание, что это mpcstate объект (не простая структура) и содержит состояния контроллера только в текущее время (не вся история до текущего времени, как структура, возвращенная sim). Он также содержит ковариационную матрицу оценщика. На каждом шаге симуляции необходимо обновлять xmpc с новыми состояниями контроллера.
xmpc % controller states
MPCSTATE object with fields
Plant: [0 0 0 0 0]
Disturbance: 0
Noise: [1x0 double]
LastMove: 0
Covariance: []
Задайте массивы рабочей области YY и UU чтобы сохранить траектории с обратной связью, чтобы они могли быть построены после симуляции.
YY=[]; UU=[];
Запустите цикл симуляции
for k=0:round(Tstop/Ts)-1 % Define measured disturbance signal v(k). You might specify a more % complex dependence on time or previous states here. v = 0; if k*Ts>=10 % raising to 1 after 10 seconds v = 1; end % Define unmeasured disturbance signal d(k). You might specify a more % complex dependence on time or previous states here. d = 0; if k*Ts>=20 % falling to -0.5 after 20 seconds d = -0.5; end % Plant equations: current output % If you have a more complex plant output behavior (including, for example, % model mismatch or nonlinearities) you might want to simulate it here. % Note that there cannot be any direct feedthrough between u and y. y = C*x + D(:,2)*v + D(:,3)*d; % calculate current output (D(:,1)=0) YY = [YY,y]; % store current output % Compute the MPC action (u) and update the controller states (xmpc) using mpcmove % Note that xmpc is updated by mpcmove even though it is an input argument! u = mpcmove(mpcobj,xmpc,y,r,v); % xmpc,y,r,v are values at current step k % Plant equations: state update % You might want to simulate a more complex plant state behavior here. x = A*x + B(:,1)*u + B(:,2)*v + B(:,3)*d; % calculate next state and update current state UU = [UU,u]; % store current input end
Постройте график результатов.
figure % create figure subplot(2,1,1) % create upper subplot axis plot(0:Ts:Tstop-Ts,YY) % plot YY versus time ylabel('y') % add a label to the upper y axis grid % add grid to upper plot title('Output') % add title to upper plot subplot(2,1,2) % create lower subplot axis plot(0:Ts:Tstop-Ts,UU) % plot UU versus time ylabel('y') % add a label to the lower y axis xlabel('Time (sec)') % add a label to the lower x axis grid % add grid to lower plot title('Input') % add title to lower plot
Если в любой момент во время симуляции вы хотите проверить оптимальные предсказанные траектории с этой точки, они могут быть возвращены mpcmove также. Предположим, что вы начинаете с текущего состояния (x и xmpc), в то время как ссылка изменяется на 0,5 и что нарушение порядка исчезло (и что оба сигнала остаются постоянными во время горизонта предсказания).
r = 0.5; % reference v = 0; % disturbance [~,info] = mpcmove(mpcobj,xmpc,y,r,v); % solve over prediction horizon
Отобразите информационную переменную.
info
info =
struct with fields:
Uopt: [11x1 double]
Yopt: [11x1 double]
Xopt: [11x6 double]
Topt: [11x1 double]
Slack: 0
Iterations: 1
QPCode: 'feasible'
Cost: 0.1399
info - структура, содержащая предсказанные оптимальные последовательности манипулируемых переменных, состояний объекта и выходов на горизонте предсказания. Эта последовательность вычисляется вместе с оптимальным перемещением путем решения квадратичной задачи оптимизации, чтобы минимизировать функцию затрат. Объект состояний и выходов в info результат применения оптимальной манипулированной переменной последовательности непосредственно к mpcobj.Model.Plant, в разомкнутом контуре. Поэтому этот процесс оптимизации без разомкнутого контура не эквивалентен симуляции замкнутого цикла, состоящего из объекта, оценщика и контроллера, используя либо sim команда или mpcmove итеративно в цикле for.
Извлеките предсказанные оптимальные траектории.
topt = info.Topt; % time yopt = info.Yopt; % predicted optimal plant model outputs uopt = info.Uopt; % predicted optimal mv sequence
Постройте график траекторий с помощью ступенчатых графиков. Ступенчатые графики используются, потому что нормальный график просто соединяет каждую вычисленную точку со следующей с помощью прямой прямой линии, в то время как оптимальные сигналы дискретного времени мгновенно скачут на каждом шаге и остаются постоянными до следующего шага. Это различие особенно заметно, поскольку для выходного сигнала объекта нужно построить только 10 интервалов (и только 3 для управляемой переменной).
figure % create new figure subplot(2,1,1) % create upper subplot stairs(topt,yopt) % plot yopt in a stairstep graph title('Optimal sequence of predicted outputs') % add title to upper subplot grid % add grid to upper subplot subplot(2,1,2) % create lower subplot stairs(topt,uopt) % plot uopt in a stairstep graph title('Optimal sequence of manipulated variables') % add title to upper subplot grid % add grid to upper subplot
Линейное представление контроллера MPC
Когда ограничения не активны, контроллер MPC ведет себя как линейный контроллер. Обратите внимание, что для конечной задачи Линейного Квадратичного Регулятора без ограничений по времени с конечным невосстановляющимся горизонтом, функция ценности зависит от времени, что заставляет оптимальный коэффициент усиления обратной связи изменяться во времени. Напротив, в MPC горизонт имеет постоянную длину, потому что он всегда отступает, приводя к инвариантной по времени функции ценности, и, следовательно, к инвариантному по времени оптимальному усилению обратной связи.
Вы можете получить форму пространства состояний контроллера MPC.
LTI = ss(mpcobj,'rv'); % get state space representation
Получите матрицы пространства состояний, чтобы симулировать линеаризированный контроллер.
[AL,BL,CL,DL] = ssdata(LTI); % get state space matrices
Инициализируйте переменные для симуляции замкнутого цикла как исходного контроллера MPC без ограничений, так и линеаризированного контроллера
mpcobj.MV = []; % remove input constraints mpcobj.OV = []; % remove output constraints Tstop = 5; % set simulation time y = 0; % set nitial measured output r = 1; % set output reference set-point (constant) u = 0; % set previous (initial) input command x = [0 0 0 0 0]'; % set initial state of plant xmpc = mpcstate(mpcobj); % set initial state of unconstrained MPC controller xL = zeros(size(BL,1),1); % set initial state of linearized MPC controller YY = []; % define workspace array to store plant outputs
Assuming no disturbance added to measured output channel #1. -->The "Model.Noise" property of the "mpc" object is empty. Assuming white noise on each measured output channel.
Симулируйте оба контроллера в замкнутом цикле с одной и той же моделью объекта управления
for k = 0:round(Tstop/Ts)-1 YY = [YY,y]; % store current output for plotting purposes % Define measured disturbance signal v(k). v = 0; if k*Ts>=10 v = 1; % raising to 1 after 10 seconds end % Define unmeasured disturbance signal d(k). d = 0; if k*Ts>=20 d = -0.5; % falling to -0.5 after 20 seconds end % Compute the control actions of both (unconstrained) MPC and linearized MPC uMPC = mpcmove(mpcobj,xmpc,y,r,v); % unconstrained MPC (this also updates xmpc) u = CL*xL+DL*[y;r;v]; % unconstrained linearized MPC % Compare the two control actions dispStr(k+1) = {sprintf('t=%5.2f, u=%7.4f (provided by LTI), u=%7.4f (provided by MPCOBJ)',k*Ts,u,uMPC)}; %#ok<*SAGROW> % Update state of unconstrained linearized MPC controller xL = AL*xL + BL*[y;r;v]; % Update plant state x = A*x + B(:,1)*u + B(:,2)*v + B(:,3)*d; % Calculate plant output y = C*x + D(:,1)*u + D(:,2)*v + D(:,3)*d; % D(:,1)=0 end
Отображение символьных массивов, содержащих действия управления
for k=0:round(Tstop/Ts)-1 disp(dispStr{k+1}); % display each string as k increases end
t= 0.00, u= 5.2478 (provided by LTI), u= 5.2478 (provided by MPCOBJ) t= 0.20, u= 3.0134 (provided by LTI), u= 3.0134 (provided by MPCOBJ) t= 0.40, u= 0.2281 (provided by LTI), u= 0.2281 (provided by MPCOBJ) t= 0.60, u=-0.9952 (provided by LTI), u=-0.9952 (provided by MPCOBJ) t= 0.80, u=-0.8749 (provided by LTI), u=-0.8749 (provided by MPCOBJ) t= 1.00, u=-0.2022 (provided by LTI), u=-0.2022 (provided by MPCOBJ) t= 1.20, u= 0.4459 (provided by LTI), u= 0.4459 (provided by MPCOBJ) t= 1.40, u= 0.8489 (provided by LTI), u= 0.8489 (provided by MPCOBJ) t= 1.60, u= 1.0192 (provided by LTI), u= 1.0192 (provided by MPCOBJ) t= 1.80, u= 1.0511 (provided by LTI), u= 1.0511 (provided by MPCOBJ) t= 2.00, u= 1.0304 (provided by LTI), u= 1.0304 (provided by MPCOBJ) t= 2.20, u= 1.0053 (provided by LTI), u= 1.0053 (provided by MPCOBJ) t= 2.40, u= 0.9920 (provided by LTI), u= 0.9920 (provided by MPCOBJ) t= 2.60, u= 0.9896 (provided by LTI), u= 0.9896 (provided by MPCOBJ) t= 2.80, u= 0.9925 (provided by LTI), u= 0.9925 (provided by MPCOBJ) t= 3.00, u= 0.9964 (provided by LTI), u= 0.9964 (provided by MPCOBJ) t= 3.20, u= 0.9990 (provided by LTI), u= 0.9990 (provided by MPCOBJ) t= 3.40, u= 1.0002 (provided by LTI), u= 1.0002 (provided by MPCOBJ) t= 3.60, u= 1.0004 (provided by LTI), u= 1.0004 (provided by MPCOBJ) t= 3.80, u= 1.0003 (provided by LTI), u= 1.0003 (provided by MPCOBJ) t= 4.00, u= 1.0001 (provided by LTI), u= 1.0001 (provided by MPCOBJ) t= 4.20, u= 1.0000 (provided by LTI), u= 1.0000 (provided by MPCOBJ) t= 4.40, u= 0.9999 (provided by LTI), u= 0.9999 (provided by MPCOBJ) t= 4.60, u= 1.0000 (provided by LTI), u= 1.0000 (provided by MPCOBJ) t= 4.80, u= 1.0000 (provided by LTI), u= 1.0000 (provided by MPCOBJ)
Постройте график результатов.
figure % create figure plot(0:Ts:Tstop-Ts,YY) % plot plant outputs grid % add grid title('Unconstrained MPC control: Plant output') % add title xlabel('Time (sec)') % add label to x axis ylabel('y') % add label to y axis
Выполнение симуляции с обратной связью, в которой все ограничения контроллера выключены, легче использовать sim, так как вам просто нужно задать ' off ' в Constraint поле связанной mpcsimopt объект опции симуляции.
SimOptions = mpcsimopt; % create simulation options object SimOptions.Constraints = 'off'; % remove all MPC constraints SimOptions.UnmeasuredDisturbance = d; % unmeasured input disturbance sim(mpcobj,Nf,r,v,SimOptions); % run closed loop simulation
Симулируйте с использованием Simulink
®Simulink является графическим блоком схемой, окружением для многодоменных системных симуляций. Можно соединить блоки, представляющие динамические системы (в этом случае объект и контроллер MPC), и моделировать замкнутый цикл. Помимо того, что система и ее взаимосвязи представлены визуально, одним из ключевых различий является то, что Simulink может использовать более широкую область значений алгоритмов интегрирования, чтобы решить базовые дифференциальные уравнения системы. Simulink также может автоматически сгенерировать код C (или PLC) и развернуть его для приложений реального времени.
% To run this part of the example, Simulink(R) is required. Check that % Simulink is installed, otherwise display a message and return if ~mpcchecktoolboxinstalled('simulink') % if Simulink not installed disp('Simulink(R) is required to run this part of the example.') return % end execution here end
Воссоздайте исходный объект mpc с исходным ограничением и исходным оценщиком по умолчанию, чтобы вы могли легко сравнить результаты.
mpcobj = mpc(plantDSS,Ts,10,3); mpcobj.MV = struct('Min',0,'Max',1,'RateMin',-10,'RateMax',10); mpcobj.Model.Disturbance = tf(sqrt(1000),[1 0]);
-->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.
Извлеките матрицы пространства состояний (непрерывного времени) объекта, которые будут управляться, поскольку они необходимы модели Simulink, которая должна быть открыта.
[A,B,C,D] = ssdata(plantCSS); % get state-space realization
Откройте ранее существовавшую модель Simulink для симуляции замкнутой системы. Модель объекта управления реализована с непрерывным блоком состояния пространство, и алгоритм интегрирования ode45 используется для вычисления его непрерывной временной характеристики.
%The block inputs are u(t), v(t) and d(t) representing the manipulated % variable, measured input disturbance, and unmeasured input disturbance, % respectively, while y(t) is the measured output. The block parameters % are the matrices forming the state space realization of the continuous-time % plant, and the initial conditions for the 5 states. % The MPC controller is implemented with an MPC controller block, which % has the workspace mpc object |mpcobj| as parameter, the manipulated % variable as output, and then the measured plant output, reference signal, % and measured plant input disturbance as inputs, respectively. % The four Scopes blocks plot the five loop signals, which are also saved % (except the reference signal) by four To-Workspace blocks. open_system('mpc_miso')
Симулируйте систему с помощью sim simulink команда. Обратите внимание, что эта команда (которая симулирует модель Simulink и эквивалентна нажатию кнопки «Запустить» в модели) отличается от sim команда, предоставляемая тулбоксом mpc (который вместо этого моделирует контроллер MPC в цикле с объектом LTI).
sim('mpc_miso')
Assuming no disturbance added to measured output channel #1. -->The "Model.Noise" property of the "mpc" object is empty. Assuming white noise on each measured output channel.
Графики в Scopes блоки эквивалентны таковым на рисунках 1,2 и 3, с незначительными различиями из-за того, что в первой симуляции (Фигуры 1 и 2) неизмеренный нарушением порядка сигнал был нулем, и то есть второй симуляции (рисунок 3) шум был добавлен к входу и выходу объекта. Также обратите внимание, что, в то время как MPC sim команда внутренне дискретизирует любую непрерывную модель объекта управления, используя ZOH, Simulink обычно использует алгоритм интегрирования (в этом примере ode45), чтобы симулировать замкнутый цикл, когда присутствует блок непрерывного времени.
Запустите симуляцию с синусоидальным выходным шумом.
Предположим, что на выходные измерения влияет синусоидальное нарушение порядка (шум одинарного датчика тона) частоты 0,1 Гц на измеренном выходе.
omega = 2*pi/10; % disturbance radial frequency
Откройте mpc_misonoise Модель Simulink. Это похоже на модель 'mpc _ miso', за исключением синусоидального нарушения порядка, вводимого на измеренном выходе. Что также отличается, так это то, что время симуляции теперь установлено на 150 секунд, и что неизмеренные и измеренные входные нарушения порядка начинают действовать через 60 и 120 секунд, соответственно. Обратите внимание, что, по-другому по сравнению с предыдущими симуляциями, неизмеренное нарушение порядка запускается первым. Это позволяет более четко построить график и проанализировать ответ на неизмеренные нарушения порядка, происходящие от 60 до 120 секунд.
open_system('mpc_misonoise') % open new Simulink model
Поскольку этот шум ожидается, установка его в модель измерения шума поможет оценщику состояния игнорировать его, тем самым улучшая качество оценок состояния и позволяя всему контроллеру лучше отклонять нарушения порядка при отслеживании выхода ссылки.
mpcobj.Model.Noise = 0.5*tf(omega^2,[1 0 omega^2]); % measurement noise model
Измените проект MPC, задав модель возмущения на неизмеренном входе в виде белого Гауссова шума с нулем среднего и отклонением 0.1. Это позволяет лучше отклонить типовые непостоянные нарушения порядка (за счет худшего отклонения входа постояннейшие нарушения порядка).
mpcobj.Model.Disturbance = tf(0.1); % static gain
Обратите внимание, что когда вы задаете статический коэффициент усиления как модель для входа нарушения порядка, выходом модель возмущения, состоящая из дискретизированного интегратора, автоматически добавляется к контроллеру. Это помогает отклонить постоянные (и медленно изменяющиеся) выходные нарушения порядка.
getoutdist(mpcobj)
-->Assuming output disturbance added to measured output channel #1 is integrated white noise.
-->A feedthrough channel in NoiseModel was inserted to prevent problems with estimator design.
ans =
A =
x1
x1 1
B =
u1
x1 0.2
C =
x1
MO1 1
D =
u1
MO1 0
Sample time: 0.2 seconds
Discrete-time state-space model.
Уменьшите вес на отслеживании выходной переменной (тем самым делая больший акцент на минимизации скорости изменения манипулируемых переменных).
mpcobj.weights = struct('MV',0,'MVRate',0.1,'OV',0.005); % new weights
Сами по себе новые веса слишком сильно штрафуют манипулируемое изменение переменной, что приводит к очень низкой пропускной способности (медленному) контроллеру. Это приведет к отставанию выходных переменных от ссылки на многие шаги. Увеличение горизонта прогнозирования до 40 позволяет контроллеру полностью принимать во внимание стоимость выходной ошибки на 40 шагах вместо только 10, тем самым делая некоторый акцент на отслеживании.
mpcobj.predictionhorizon = 40; % new prediction horizon
Запустите симуляцию в течение 145 секунд
sim('mpc_misonoise',145) % the second argument is the simulation duration
-->Assuming output disturbance added to measured output channel #1 is integrated white noise. -->A feedthrough channel in NoiseModel was inserted to prevent problems with estimator design.
Фильтр Калмана успешно учится игнорировать шум измерения через 50 секунд. Неизмеренные и измеренные нарушения порядка отклоняются в интервале времени от 10 до 20 секунд. Наконец, как и ожидалось, манипулируемая переменная остается в интервале между 0 и 1.
bdclose('all') % close all open Simulink models without saving any change close all % close all open figures
См. также
mpc | MPC Controller | MPC Designer