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

Откройте модель 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 function, убедитесь, что модель является минимальной реализацией.

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)