Открытие модели

Откройте модель Simulink.

mdl = 'Boiler_Demo';
open_system(mdl)

При открытии модели управления котлом программа инициализирует размеры контроллера. u0 и y0 устанавливаются после вычисления рабочей точки и поэтому первоначально устанавливаются равными нулю. Наблюдатель и регулятор вычисляются на этапе проектирования контроллера и также первоначально устанавливаются равными нулю.

Найти номинальную рабочую точку и модель линеаризации

Значения начального состояния модели определяются в модели Simulink. Используя эти значения состояния, найдите рабочую точку установившегося состояния, используя findop функция.

Создайте спецификацию рабочей точки, в которой известны значения состояния.

opspec = operspec(mdl);
opspec.States(1).Known = 1;
opspec.States(2).Known = 1;
opspec.States(3).Known = [1;1];

Откорректируйте спецификацию рабочей точки, чтобы указать, что входные данные должны быть вычислены и что они имеют нижний предел.

opspec.Inputs(1).Known = [0;0];     % Inputs unknown
opspec.Inputs(1).Min = [0;0];       % Input minimum value

Добавить выходную спецификацию в спецификацию рабочей точки; это необходимо для обеспечения вычисления рабочей точки вывода во время процесса решения.

opspec = addoutputspec(opspec,[mdl '/Boiler'],1);
opspec.Outputs(1).Known = 0;    % Outputs unknown
opspec.Outputs(1).Min = 0;      % Output minimum value

Вычислите операционную точку и создайте отчет о поиске операционной точки.

[opSS,opReport] = findop(mdl,opspec);

 Operating point search report:
---------------------------------

 Operating point search report for the Model Boiler_Demo.
 (Time-Varying Components Evaluated at time t=0)

Operating point specifications were successfully met.
States: 
----------
(1.) Boiler_Demo/Boiler/Steam volume
      x:           5.6      dx:      7.85e-13 (0)
(2.) Boiler_Demo/Boiler/Temperature
      x:           180      dx:     -5.93e-14 (0)
(3.) Boiler_Demo/Observer/Internal
      x:             0      dx:             0 (0)
      x:             0      dx:             0 (0)

Inputs: 
----------
(1.) Boiler_Demo/Input
      u:      2.41e+05    [0 Inf]
      u:           100    [0 Inf]

Outputs: 
----------
(1.) Boiler_Demo/Boiler
      y:         1e+03    [0 Inf]

Перед линеаризацией модели вокруг этой точки укажите входной и выходной сигналы для линейной модели. Сначала задайте входные точки для линеаризации.

Boiler_io(1) = linio([mdl '/Sum'],1,'input');
Boiler_io(2) = linio([mdl '/Demux'],2,'input');

Затем задайте выходные точки с разомкнутым контуром для линеаризации.

Boiler_io(3) = linio([mdl '/Boiler'],1,'openoutput');
setlinio(mdl,Boiler_io);

Найдите линейную модель вокруг выбранной рабочей точки.

Lin_Boiler = linearize(mdl,opSS,Boiler_io);

Наконец, использование minreal убедитесь, что модель является минимальной реализацией.

Lin_Boiler = minreal(Lin_Boiler);

1 state removed.

Проектный регулятор и государственный наблюдатель

Используя эту линейную модель, спроектируйте регулятор LQR и наблюдателя состояния фильтра Калмана. Сначала найдите смещения контроллера, чтобы убедиться, что контроллер работает вокруг выбранной точки линеаризации, путем извлечения вычисленной рабочей точки.

u0 = opReport.Inputs.u;
y0 = opReport.Outputs.y;

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

Q = diag(1e8);                  % Output regulation
R = diag([1e2,1e6]);            % Input limitation
[K,S,E] = lqry(Lin_Boiler,Q,R);

Спроектируйте государственного наблюдателя Калмана, используя kalman функция. Для этого примера основным источником шума является технологический шум. Шум поступает в систему только через один вход, следовательно, в виде G и H.

[A,B,C,D] = ssdata(Lin_Boiler);
G = B(:,1);
H = 0;
QN = 1e4;
RN = 1e-1;
NN = 0;
[Kobsv,L,P] = kalman(ss(A,[B G],C,[D H]),QN,RN);

Моделирование модели

Смоделировать модель для разработанного контроллера.

sim(mdl)

Постройте график входных и выходных сигналов процесса. На следующем рисунке показан сигнал срабатывания питательной воды в кг/с.

plot(FeedWater.time/60,FeedWater.signals.values)
title('Feedwater flow rate [kg/s]');
ylabel('Flow [kg/s]')
xlabel('time [min]')
grid on

На следующем графике показан сигнал включения тепла в кДж.

plot(Heat.time/60,Heat.signals.values/1000)
title('Applied heat [kJ]');
ylabel('Heat [kJ]')
xlabel('time [min]')
grid on

На следующем рисунке показано тепловое возмущение в кДж. Возмущение изменяется на 50% от номинального значения тепла.

plot(HeatDist.time/60,HeatDist.signals.values/1000)
title('Heat disturbance [kJ]');
ylabel('Heat [kJ]')
xlabel('time [min]')
grid on

На рисунке ниже показано соответствующее давление в емкости в кПа. Давление изменяется примерно на 1% от номинального значения, даже если возмущение относительно велико.

plot(DrumPressure.time/60,DrumPressure.signals.values)
title('Drum pressure [kPa]');
ylabel('Pressure [kPa]')
xlabel('time [min]')
grid on

bdclose(mdl)