Создайте файл MEX, который решает (типовую или многоступенчатую) нелинейную задачу MPC управления
mexFcn = buildMEX(nlobj,mexName,coreData,onlineData)nlmpcmove. Файл MEX создается в текущей рабочей папке.
Создайте нелинейный контроллер MPC с четырьмя состояниями, двумя выходами и одним входом.
nlobj = nlmpc(4,2,1);
In standard cost function, zero weights are applied by default to one or more OVs because there are fewer MVs than OVs.
Укажите шаг расчета и горизонты контроллера.
Ts = 0.1; nlobj.Ts = Ts; nlobj.PredictionHorizon = 10; nlobj.ControlHorizon = 5;
Задайте функцию состояния для контроллера, которая находится в файле pendulumDT0.m. Эта модель в дискретном времени интегрирует модель в непрерывном времени, заданную в pendulumCT0.m использование многостадийного прямого метода Эйлера.
nlobj.Model.StateFcn = "pendulumDT0";
nlobj.Model.IsContinuousTime = false;Модель предсказания использует необязательный параметр Ts для представления шага расчета. Задайте количество параметров и создайте вектор параметра.
nlobj.Model.NumberOfParameters = 1;
params = {Ts};Задайте выходную функцию модели, передав параметр времени расчета как входной параметр.
nlobj.Model.OutputFcn = "pendulumOutputFcn";Задайте стандартные ограничения для контроллера.
nlobj.Weights.OutputVariables = [3 3]; nlobj.Weights.ManipulatedVariablesRate = 0.1; nlobj.OV(1).Min = -10; nlobj.OV(1).Max = 10; nlobj.MV.Min = -100; nlobj.MV.Max = 100;
Проверьте функции модели предсказания.
x0 = [0.1;0.2;-pi/2;0.3]; u0 = 0.4; validateFcns(nlobj,x0,u0,[],params);
Model.StateFcn is OK. Model.OutputFcn is OK. Analysis of user-provided model, cost, and constraint functions complete.
Только два из состояний объекта можно измерить. Поэтому создайте расширенный фильтр Калмана для оценки четырех состояний объекта. Его функция перехода к состоянию задана в pendulumStateFcn.m и его функция измерения определена в pendulumMeasurementFcn.m.
EKF = extendedKalmanFilter(@pendulumStateFcn,@pendulumMeasurementFcn);
Задайте начальные условия для симуляции, инициализируйте расширенное состояние фильтра Калмана и задайте нуль начальный манипулированное значение переменных.
x0 = [0;0;-pi;0]; y0 = [x0(1);x0(3)]; EKF.State = x0; mv0 = 0;
Создайте структуры данных генерации кода для контроллера, определяя начальные условия и параметры.
[coreData,onlineData] = getCodeGenerationData(nlobj,x0,mv0,params);
Задайте значение выхода ссылки в интерактивной структуре данных.
onlineData.ref = [0 0];
Создайте MEX-функцию для решения нелинейной задачи MPC управления. Эта MEX-функция создаётся в текущую рабочую директорию.
mexFcn = buildMEX(nlobj,"myController",coreData,onlineData);Generating MEX function "myController" from nonlinear MPC to speed up simulation. Code generation successful. MEX function "myController" successfully generated.
Запустите симуляцию для 10 секунд. Во время каждого контрольного интервала:
Исправьте предыдущее предсказание с помощью измерения тока.
Вычислите оптимальное управление перемещается с помощью MEX-функции. Эта функция возвращает вычисленные оптимальные последовательности в onlineData. Передача обновленной структуры данных в MEX-функцию в следующем контрольном интервале предоставляет начальные предположения для оптимальных последовательностей.
Спрогнозируйте состояния модели.
Примените первое вычисленное оптимальное перемещение управления к объекту, обновляя состояния объекта.
Сгенерируйте данные датчика с белым шумом.
Сохраните состояния объекта.
mv = mv0; y = y0; x = x0; Duration = 10; xHistory = x0; for ct = 1:(Duration/Ts) % Correct previous prediction xk = correct(EKF,y); % Compute optimal control move [mv,onlineData] = myController(xk,mv,onlineData); % Predict prediction model states for the next iteration predict(EKF,[mv; Ts]); % Implement first optimal control move x = pendulumDT0(x,mv,Ts); % Generate sensor data y = x([1 3]) + randn(2,1)*0.01; % Save plant states xHistory = [xHistory x]; end
Постройте график получившихся траекторий состояния.
figure subplot(2,2,1) plot(0:Ts:Duration,xHistory(1,:)) xlabel('time') ylabel('z') title('cart position') subplot(2,2,2) plot(0:Ts:Duration,xHistory(2,:)) xlabel('time') ylabel('zdot') title('cart velocity') subplot(2,2,3) plot(0:Ts:Duration,xHistory(3,:)) xlabel('time') ylabel('theta') title('pendulum angle') subplot(2,2,4) plot(0:Ts:Duration,xHistory(4,:)) xlabel('time') ylabel('thetadot') title('pendulum velocity')
getCodeGenerationData | nlmpc | nlmpcmove | nlmpcmoveCodeGeneration