В этом примере показано, как сконструировать прогнозирующий контроллер модели для установки с одним измеренным выходом, одной управляемой переменной, одним измеренным возмущением и одним неизмеренным возмущением.
Определите модель завода. В этом примере используются функции непрерывного переноса времени с каждого входа на выход. Отображение модели завода в командной строке.
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 автоматически преобразует установку непрерывного времени (в любом формате) в модель пространства состояния дискретного времени, используя Удержание нулевого порядка (Zero 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
Создайте объект контроллера, указав время выборки, а также горизонты прогнозирования и управления (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 Команда 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 шагов с отсчетом r и измеренным входным возмущением 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, для которой невозможно рассчитать СН.
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 горизонт имеет постоянную длину, потому что он всегда отступает, что приводит к функции инвариантного значения времени и, следовательно, к оптимальному коэффициенту усиления обратной связи.
Вы можете получить форму «state-space» контроллера 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 может использовать более широкий диапазон алгоритмов интеграции для решения лежащих в основе системы дифференциальных уравнений. 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 и эквивалентная нажатию кнопки «Выполнить» в модели) отличается от команды 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 | Разработчик MPC