Руководство: линейное проектирование

Руководство: линейное проектирование

Линейная модель

В этом примере в качестве объекта взят смесительный реактор непрерывного действия (CSTR) с впусками холодной и горячей воды. Воды вытекает из отверстия на дне бака. Обрабатываемыми входами являются расход холодной ( ) и горячей ( ) воды, а измеряемыми выходами — уровень ( ) и температура ( ) жидкости:

На следующем рисунке изображены устройства, установленные на реакторе CSTR:

В рабочих точках с установившимся состоянием:

следующая линейная модель точно описывает динамику объекта:

Сначала с помощью функции setop! необходимо создать объект LinModel для операций с рабочими точками:

using ModelPredictiveControl, ControlSystemsBase
G = [ tf(1.90, [18, 1]) tf(1.90, [18, 1]);
      tf(-0.74,[8, 1])  tf(0.74, [8, 1]) ]
Ts = 2.0
model = setop!(LinModel(G, Ts), uop=[20, 20], yop=[50, 30])

Discrete-time linear model with a sample time Ts = 2.0 s and:
 2 manipulated inputs u
 2 states x
 2 outputs y
 0 measured disturbances d

Объект model будет служить двум целям: построение контроллера и моделирование объекта для тестирования проекта. Время дискретизации составляет 2 секунды, так что период управления будет также иметь длительность 2 секунды.

Прогнозирующий контроллер линейной модели

Прогнозирующий контроллер линейной модели (MPC) будет контролировать как уровень воды , так и температуру в баке. Кроме того, уровень в баке не должен падать ниже 45:

Для проектирования контроллеров LinMPC мы включаем линейное ограничение уровня с помощью функции setconstraint! (если ограничения нет, следует использовать значения ±Inf):

mpc = LinMPC(model, Hp=10, Hc=2, Mwt=[1, 1], Nwt=[0.1, 0.1])
mpc = setconstraint!(mpc, ymin=[45, -Inf])

LinMPC controller with a sample time Ts = 2.0 s, OSQP optimizer, SteadyKalmanFilter estimator and:
 10 prediction steps Hp
  2 control steps Hc
  2 manipulated inputs u (0 integrating states)
  4 estimated states x̂
  2 measured outputs ym (2 integrating states)
  0 unmeasured outputs yu
  0 measured disturbances d

Здесь именованные аргументы Hp и Hc — это соответственно горизонты прогнозирования и управления, а Mwt и Nwt — веса отслеживания уставки выхода и гашения перемещений. По умолчанию контроллеры LinMPC используют OSQP для решения задачи, мягкие ограничения для прогнозов на выходах с целью обеспечения применимости и SteadyKalmanFilter для оценки состояний объекта [^1]. Внимательный читатель мог также обратить внимание, что фильтр Калмана оценивает два дополнительных состояния по сравнению с моделью объекта. Это интегрирующие состояния для неизмеряемых возмущений объекта, которые по умолчанию автоматически добавляются к выходам модели, если сохраняется наблюдаемость (в описании типа SteadyKalmanFilter приводятся подробные сведения). [^1]: В качестве альтернативы наблюдателю состояния мы могли бы использовать структуру InternalModel с mpc = LinMPC(InternalModel(model), Hp=15, Hc=2, Mwt=[1, 1], Nwt=[0.1, 0.1]). Она была протестирована на примере, приведенном на этой странице, и показала схожие результаты.

Перед замыканием контура мы вызываем функцию initstate! с фактическими входами и измеренными значениями объекта, чтобы обеспечить мягкую передачу управления. Так как объект моделируется с помощью model, его выход инициализирует состояния. С этой целью объекты LinModel являются вызываемыми (псевдоним для evaloutput):

u, y = model.uop, model() # или, что то же самое: y = evaloutput(model)
initstate!(mpc, u, y)

После этого можно замкнуть контур и протестировать производительность mpc в симуляторе, применив шаговые изменения к уставкам выходов и возмущению по нагрузке :

function test_mpc(mpc, model)
    N = 200
    ry, ul = [50, 30], 0
    u_data, y_data, ry_data = zeros(model.nu, N), zeros(model.ny, N), zeros(model.ny, N)
    for i = 1:N
        i == 51  && (ry = [50, 35])
        i == 101 && (ry = [54, 30])
        i == 151 && (ul = -20)
        y = model() # моделируемые измерения
        u = mpc(ry) # или, что то же самое: u = moveinput!(mpc, ry)
        u_data[:,i], y_data[:,i], ry_data[:,i] = u, y, ry
        updatestate!(mpc, u, y) # обновляем оценку состояния mpc
        updatestate!(model, u + [0; ul]) # обновляем симулятор, внося возмущение по нагрузке
    end
    return u_data, y_data, ry_data
end
u_data, y_data, ry_data = test_mpc(mpc, model)
t_data = Ts*(0:(size(y_data,2)-1))

Объекты LinMPC также могут вызываться в качестве альтернативы синтаксису moveinput!. Стоит отметить, что дополнительные сведения, такие как оптимальные прогнозы на выходах , можно получить путем вызова getinfo после решения задачи. Кроме того, при вызове updatestate! для объекта mpc обновляется его внутреннее состояние для следующего периода управления (это сделано преднамеренно; в разделе Функции: средства оценки состояний приводится обоснование). Вот почему этот вызов совершается в конце цикла for. Та же логика действительна для model.

Наконец, мы строим график результатов теста с замкнутым контуром с помощью пакета Plots:

using Plots
function plot_data(t_data, u_data, y_data, ry_data)
    p1 = plot(t_data, y_data[1,:], label="meas."); ylabel!("level")
    plot!(p1, t_data, ry_data[1,:], label="setpoint", linestyle=:dash, linetype=:steppost)
    plot!(p1, t_data, fill(45,size(t_data)), label="min", linestyle=:dot, linewidth=1.5)
    p2 = plot(t_data, y_data[2,:], label="meas.", legend=:topleft); ylabel!("temp.")
    plot!(p2, t_data, ry_data[2,:],label="setpoint", linestyle=:dash, linetype=:steppost)
    p3 = plot(t_data,u_data[1,:],label="cold", linetype=:steppost); ylabel!("flow rate")
    plot!(p3, t_data,u_data[2,:],label="hot", linetype=:steppost); xlabel!("time (s)")
    return plot(p1, p2, p3, layout=(3,1))
end
plot_data(t_data, u_data, y_data, ry_data)

По сравнению с конфигурацией по умолчанию добавление интегрирующих состояний на входах модели может повысить производительность при замкнутом контуре. Возмущения по нагрузке на самом деле очень часто встречаются в реальных задачах управления. Создание LinMPC с интеграторами входов:

mpc2 = LinMPC(model, Hp=10, Hc=2, Mwt=[1, 1], Nwt=[0.1, 0.1], nint_u=[1, 1])
mpc2 = setconstraint!(mpc2, ymin=[45, -Inf])

LinMPC controller with a sample time Ts = 2.0 s, OSQP optimizer, SteadyKalmanFilter estimator and:
 10 prediction steps Hp
  2 control steps Hc
  2 manipulated inputs u (2 integrating states)
  4 estimated states x̂
  2 measured outputs ym (0 integrating states)
  0 unmeasured outputs yu
  0 measured disturbances d

позволяет ускорить исключение возмущения по нагрузке и предотвратить нарушение ограничения уровня:

setstate!(model, zeros(model.nx))
u, y = model.uop, model()
initstate!(mpc2, u, y)
u_data, y_data, ry_data = test_mpc(mpc2, model)
plot_data(t_data, u_data, y_data, ry_data)

Добавление компенсации по каналу прямой связи

Предположим, что возмущение по нагрузке в последней части на самом деле вызывается отдельным трубопроводом, по которому горячая вода поступает в бак. Измерение этого расхода позволяет включить в контроллер компенсацию по каналу прямой связи. Новая модель объекта выглядит так:

Нам необходимо создать новую модель LinModel, которая включает измеряемое возмущение и рабочую точку :

model_d = setop!(LinModel([G G[1:2, 2]], Ts, i_d=[3]), uop=[20, 20], yop=[50, 30], dop=[20])

Discrete-time linear model with a sample time Ts = 2.0 s and:
 2 manipulated inputs u
 4 states x
 2 outputs y
 1 measured disturbances d

На основе этой модели создается контроллер LinMPC:

mpc_d = LinMPC(model_d, Hp=10, Hc=2, Mwt=[1, 1], Nwt=[0.1, 0.1])
mpc_d = setconstraint!(mpc_d, ymin=[45, -Inf])

LinMPC controller with a sample time Ts = 2.0 s, OSQP optimizer, SteadyKalmanFilter estimator and:
 10 prediction steps Hp
  2 control steps Hc
  2 manipulated inputs u (0 integrating states)
  6 estimated states x̂
  2 measured outputs ym (2 integrating states)
  0 unmeasured outputs yu
  1 measured disturbances d

Кроме того, требуется новая тестовая функция, которая передает измеряемое возмущение в контроллер:

function test_mpc_d(mpc_d, model)
    N = 200
    ry, ul = [50, 30], 0
    dop = mpc_d.estim.model.dop
    u_data, y_data, ry_data = zeros(model.nu, N), zeros(model.ny, N), zeros(model.ny, N)
    for i = 1:N
        i == 51  && (ry = [50, 35])
        i == 101 && (ry = [54, 30])
        i == 151 && (ul = -20)
        d = ul .+ dop   # моделируемое измеряемое возмущение
        y = model()     # моделируемые измерения
        u = mpc_d(ry, d) # кроме того, передаем в контроллер измеряемое возмущение d
        u_data[:,i], y_data[:,i], ry_data[:,i] = u, y, ry
        updatestate!(mpc_d, u, y, d)    # обновляем оценку после введения измеряемого возмущения d
        updatestate!(model, u + [0; ul]) # обновляем симулятор
    end
    return u_data, y_data, ry_data
end

Новая компенсация по каналу прямой связи позволяет почти идеально исключить возмущение по нагрузке:

setstate!(model, zeros(model.nx))
u, y, d = model.uop, model(), mpc_d.estim.model.dop
initstate!(mpc_d, u, y, d)
u_data, y_data, ry_data = test_mpc_d(mpc_d, model)
plot_data(t_data, u_data, y_data, ry_data)

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