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

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

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

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

Система представлена на следующем рисунке:

Модель объекта нелинейная:

Здесь — это ускорение свободного падения в м/с², — длина маятника в м, — коэффициент трения в точке вращения в кг/с, а — масса груза на конце маятника в кг. Система дискретизируется с помощью явного метода Эйлера с целью создания модели NonLinModel:

using ModelPredictiveControl
function pendulum(par, x, u)
    g, L, K, m = par        # (м/с²), (м), (кг/с), (кг)
    θ, ω = x[1], x[2]       # (рад), (рад/с)
    τ  = u[1]               # (Н м)
    dθ = ω
    dω = -g/L*sin(θ) - K/m*ω + τ/m/L^2
    return [dθ, dω]
end
# объявляем константы для скорости, чтобы избежать нестабильности типов в функции f:
const par, Ts = (9.8, 0.4, 1.2, 0.3), 0.1
f(x, u, _ ) = x + Ts*pendulum(par, x, u) # Метод Эйлера
h(x, _ )    = [180/π*x[1]]  # [°]
nu, nx, ny = 1, 2, 1
model = NonLinModel(f, h, Ts, nu, nx, ny)

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

Выходная функция преобразует угол в градусы. Обратите внимание, что для ввода специальных символов, таких как , в REPL Julia или VS Code, можно ввести \theta и нажать клавишу <TAB>. Кортеж par и переменная Ts объявлены здесь как константы для лучшей производительности. Чтобы быстро проверить работоспособность, рекомендуется сначала смоделировать model с помощью sim!:

using Plots
u = [0.5]
N = 35
res = sim!(model, N, u)
plot(res, plotu=false)

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

Объект UnscentedKalmanFilter оценивает состояние объекта:

σQ=[0.1, 0.5]; σR=[0.5]; nint_u=[1]; σQint_u=[0.1]
estim = UnscentedKalmanFilter(model; σQ, σR, nint_u, σQint_u)

UnscentedKalmanFilter estimator with a sample time Ts = 0.1 s, NonLinModel and:
 1 manipulated inputs u (1 integrating states)
 3 estimated states x̂
 1 measured outputs ym (0 integrating states)
 0 unmeasured outputs yu
 0 measured disturbances d

Векторы σQ и σR — это среднеквадратичные отклонения шумов процесса и датчика соответственно. Значение скорости здесь больше (второе значение σQ), так как в уравнении есть неопределенный параметр: коэффициент трения . Кроме того, аргумент nint_u явным образом добавляет на входе модели одно интегрирующее состояние, крутящий момент двигателя , со связанным среднеквадратичным отклонением σQint_u, равным 0,1 Н м. Настройка оценки проверяется на объекте с коэффициентом трения больше на 25 %:

const par_plant = (par[1], par[2], 1.25*par[3], par[4])
f_plant(x, u, _) = x + Ts*pendulum(par_plant, x, u)
plant = NonLinModel(f_plant, h, Ts, nu, nx, ny)
res = sim!(estim, N, [0.5], plant=plant, y_noise=[0.5])
plot(res, plotu=false, plotxwithx̂=true)

Оценка — это интегрирующее состояние по крутящему моменту , которое компенсирует статические погрешности. Производительность фильтра Калмана кажется достаточной для управления.

Так как крутящий момент двигателя ограничен диапазоном --1,5—1,5 Н м, в NonLinMPC включаются входные ограничения:

nmpc = NonLinMPC(estim, Hp=20, Hc=2, Mwt=[0.5], Nwt=[2.5])
nmpc = setconstraint!(nmpc, umin=[-1.5], umax=[+1.5])

NonLinMPC controller with a sample time Ts = 0.1 s, Ipopt optimizer, UnscentedKalmanFilter estimator and:
 20 prediction steps Hp
  2 control steps Hc
  1 manipulated inputs u (1 integrating states)
  3 estimated states x̂
  1 measured outputs ym (0 integrating states)
  0 unmeasured outputs yu
  0 measured disturbances d

Для тестирования производительности mpc применительно к объекту plant мы применяем угловую уставку 180° (обратное положение):

res_ry = sim!(nmpc, N, [180.0], plant=plant, x0=[0, 0], x̂0=[0, 0, 0])
plot(res_ry)

Контроллер кажется достаточно устойчивым к колебаниям коэффициента . Начиная с этого обратного положения реакция замкнутого контура на скачкообразные возмущения в 10° также удовлетворительная:

res_yd = sim!(nmpc, N, [180.0], plant=plant, x0=[π, 0], x̂0=[π, 0, 0], y_step=[10])
plot(res_yd)

Экономичный прогнозирующий контроллер модели

Экономичный контроллер MPC служит той же цели, но позволяет снизить затраты. В данном практическом случае контроллер будет пытаться снизить энергопотребление двигателя. Мощность (W), передаваемая двигателем на маятник, описывается формулой:

Таким образом, работа (J), выполняемая двигателем за время с по , описывается так:

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

Целевая функция теперь будет содержать аддитивный член, который накладывает штраф на работу, выполняемую электродвигателем , с целью снижения энергопотребления. Обратите внимание, что — это функция от обрабатываемого входа и угловой скорости , состояния, которое не измеряется (измеряется здесь только угол ). Так как аргументы экономичной функции NonLinMPC JE не включают в себя состояния, скорость теперь определяется как неизмеряемый выход для синтезирования фильтра Калмана так же, как в предыдущем случае ( и ):

h2(x, _ ) = [180/π*x[1], x[2]]
nu, nx, ny = 1, 2, 2
model2 = NonLinModel(f      , h2, Ts, nu, nx, ny)
plant2 = NonLinModel(f_plant, h2, Ts, nu, nx, ny)
estim2 = UnscentedKalmanFilter(model2; σQ, σR, nint_u, σQint_u, i_ym=[1])

UnscentedKalmanFilter estimator with a sample time Ts = 0.1 s, NonLinModel and:
 1 manipulated inputs u (1 integrating states)
 3 estimated states x̂
 1 measured outputs ym (0 integrating states)
 1 unmeasured outputs yu
 0 measured disturbances d

Кроме того, требуется объект plant2 на основе h2, так как функция sim! ожидает, что выходной вектор в аргументе plant соответствует выходному вектору модели в аргументе mpc. Теперь можно определить функцию и контроллер empc:

function JE(UE, ŶE, _ )
    τ, ω = UE[1:end-1], ŶE[2:2:end-1]
    return Ts*sum(τ.*ω)
end
empc = NonLinMPC(estim2, Hp=20, Hc=2, Mwt=[0.5, 0], Nwt=[2.5], Ewt=4.5e3, JE=JE)
empc = setconstraint!(empc, umin=[-1.5], umax=[+1.5])

NonLinMPC controller with a sample time Ts = 0.1 s, Ipopt optimizer, UnscentedKalmanFilter estimator and:
 20 prediction steps Hp
  2 control steps Hc
  1 manipulated inputs u (1 integrating states)
  3 estimated states x̂
  1 measured outputs ym (0 integrating states)
  1 unmeasured outputs yu
  0 measured disturbances d

Именованный аргумент Ewt взвешивает затраты относительно других членов в целевой функции. Член должен иметь достаточно большое значение, чтобы иметь значимость, но слишком большое значение может привести к статической погрешности угловой уставки. Второй элемент Mwt равен нулю, так как для отслеживания уставки скорость не запрашивается. Реакция замкнутого контура на уставку 180° аналогичная:

res2_ry = sim!(empc, N, [180, 0], plant=plant2, x0=[0, 0], x̂0=[0, 0, 0])
plot(res2_ry)

а энергопотребление немного ниже:

function calcW(res)
    τ, ω = res.U_data[1, 1:end-1], res.X_data[2, 1:end-1]
    return Ts*sum(τ.*ω)
end
Dict(:W_nmpc => calcW(res_ry), :W_empc => calcW(res2_ry))

Dict{Symbol, Float64} with 2 entries:
  :W_empc => 4.19004
  :W_nmpc => 4.22099

Кроме того, при ступенчатом возмущении 10°:

res2_yd = sim!(empc, N, [180; 0]; plant=plant2, x0=[π, 0], x̂0=[π, 0, 0], y_step=[10, 0])
plot(res2_yd)

новый контроллер может рекуперировать больше энергии маятника (то есть выполнять отрицательную работу):

Dict(:W_nmpc => calcW(res_yd), :W_empc => calcW(res2_yd))

Dict{Symbol, Float64} with 2 entries:
  :W_empc => -0.13526
  :W_nmpc => -0.0844623

Естественно, это даст выигрыш только в том случае, если в электронном блоке двигателя есть цепь, отвечающая за рекуперацию.

Линеаризация модели

Нелинейный контроллер MPC требует больших вычислительных затрат, чем LinMPC. В реальных условиях задача всегда должна решаться за время меньшее, чем время дискретизации с. В электронных или механических системах это требование иногда трудно выполнить из-за быстрого протекания процессов. Чтобы упростить проектирование и сравнение с LinMPC, функция linearize поддерживает автоматическую линеаризацию модели NonLinModel на основе ForwardDiff.jl. Сначала мы линеаризуем model в точке рад и (обратное положение):

linmodel = linearize(model, x=[π, 0], u=[0])

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

Стоит отметить, что метод Эйлера не лучшим образом подходит для линеаризации объекта model из-за низкой точности (аппроксимация плохой аппроксимации). Фильтр SteadyKalmanFilter и контроллер LinMPC создаются на основе модели linmodel:

kf  = SteadyKalmanFilter(linmodel; σQ, σR, nint_u, σQint_u)
mpc = LinMPC(kf, Hp=20, Hc=2, Mwt=[0.5], Nwt=[2.5])
mpc = setconstraint!(mpc, umin=[-1.5], umax=[+1.5])

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

Линейный контроллер испытывает трудности с исключением ступенчатого возмущения 10°:

res_lin = sim!(mpc, N, [180.0]; plant, x0=[π, 0], y_step=[10])
plot(res_lin)

В такой ситуации может помочь решение задачи оптимизации mpc с помощью оптимизатора DAQP, а не решателя по умолчанию OSQP. В документации указано, что DAQP может иметь более высокую производительность при работе с матрицами небольшой или средней плотности и неустойчивыми полюсами ^[1], как в данном случае (абсолютное значение неустойчивых полюсов больше единицы):

embedded quadratic programming using recursive LDLᵀ updates. IEEE Trans. Autom. Contr., 67(8). https://doi.org/doi:10.1109/TAC.2022.3176430.

using LinearAlgebra; poles = eigvals(linmodel.A)

2-element Vector{Float64}:
 0.26614608739843437
 1.3338539126015654

Для установки решателя выполните следующую команду:

using Pkg; Pkg.add("DAQP")

Создание LinMPC с оптимизатором DAQP:

using JuMP, DAQP
daqp = Model(DAQP.Optimizer, add_bridges=false)
mpc2 = LinMPC(kf, Hp=20, Hc=2, Mwt=[0.5], Nwt=[2.5], optim=daqp)
mpc2 = setconstraint!(mpc2, umin=[-1.5], umax=[+1.5])

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

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

res_lin2 = sim!(mpc2, N, [180.0]; plant, x0=[π, 0], y_step=[10])
plot(res_lin2)

Производительность замкнутого контура по-прежнему ожидаемо ниже, чем у нелинейного контроллера, но вычисления производятся примерно в 2000 раз быстрее (в среднем 0,00002 с на временной шаг вместо 0,04 с). Обратите внимание, что модель linmodel действительна только для угловых положений, близких к 180°. При больших отклонениях от этой рабочей точки требуется несколько линеаризованных моделей и контроллеров. Это называется табличным управлением усилением.

1. Arnström, D., Bemporad, A., and Axehill, D. (2022). A dual active-set solver for