Этот пример показывает, как симулировать и сгенерировать код для прогнозирующего контроллера модели, который использует пользовательский квадратичный решатель программирования (QP). Объект для этого примера является сервоприводом постоянного тока в Simulink ®.
Модель электродвигателя является линейной динамической системой, описанной в [1]. plant - модель пространства состояний двигателя в непрерывном времени. tau - максимально допустимый крутящий момент, который вы используете в качестве выхода ограничения.
[plant,tau] = mpcmotormodel;
Объект имеет один вход, входное напряжение двигателя. Контроллер MPC использует этот вход как манипулируемую переменную (MV). Объект имеет два выхода, угловое положение двигателя и крутящий момент вала. Угловое положение является измеренным выходом (MO), и крутящий момент вала не измеряется (UO).
plant = setmpcsignals(plant,'MV',1,'MO',1,'UO',2);
Ограничьте управляемую переменную, чтобы она находилась между +/- 220 вольт. Поскольку входы и выходные параметры объекта имеют различные порядки величины, для облегчения настройки используйте шкалу коэффициенты. Типичными вариантами шкалы коэффициента являются верхний/нижний предел или рабочая область значений.
MV = struct('Min',-220,'Max',220,'ScaleFactor',440);
Нет ограничений на угловое положение. Задайте верхнюю и нижнюю границы крутящего момента на валу во время первых трех шагов предсказания горизонта. Чтобы определить эти границы, используйте tau.
OV = struct('Min',{-Inf, [-tau;-tau;-tau;-Inf]},... 'Max',{Inf, [tau;tau;tau;Inf]},'ScaleFactor',{2*pi, 2*tau});
Задача управления состоит в том, чтобы достичь нулевой ошибки отслеживания для углового положения. Поскольку у вас есть только одна управляемая переменная, позволите крутящему моменту вала поплавать в пределах своего ограничения, установив его настройочный вес равным нулю.
Weights = struct('MV',0,'MVRate',0.1,'OV',[0.1 0]);
Укажите шаг расчета и горизонты и создайте контроллер MPC, используя plant как прогнозирующая модель.
Ts = 0.1; % Sample time p = 10; % Prediction horizon m = 2; % Control horizon mpcobj = mpc(plant,Ts,p,m,Weights,MV,OV);
Чтобы запустить оставшийся пример, требуется Simulink.
if ~mpcchecktoolboxinstalled('simulink') disp('Simulink is required to run this example.') return end
Откройте модель Simulink, которая имитирует управление с обратной связью двигателя постоянного сервопривода с помощью контроллера MPC. По умолчанию MPC использует встроенный решатель QP, который использует алгоритм KWIK.
mdl = 'mpc_customQPcodegen';
open_system(mdl)
Запустите симуляцию
sim(mdl)
-->Converting model to discrete time. Assuming no disturbance added to measured output channel #1. -->The "Model.Noise" property of the "mpc" object is empty. Assuming white noise on each measured output channel.
Сохраните входные и выходные сигналы объекта в рабочем пространстве MATLAB.
uKWIK = u; yKWIK = y;
Чтобы изучить, как пользовательский решатель ведет себя в тех же условиях, включите пользовательский решатель в контроллере MPC.
mpcobj.Optimizer.CustomSolver = true;
Вы также должны предоставить функцию MATLAB ®, которая удовлетворяет следующим требованиям:
Имя функции должно быть mpcCustomSolver.
Входные и выходные аргументы должны совпадать с аргументами в файле шаблона.
Функция должна находиться в пути MATLAB.
В этом примере используйте пользовательский решатель QP, заданный в файле шаблона mpcCustomSolverCodeGen_TemplateEML.txt, который реализует dantzig алгоритм и подходит для генерации кода. Сохраните функцию в рабочей папке следующим mpcCustomSolver.m.
src = which('mpcCustomSolverCodeGen_TemplateEML.txt'); dest = fullfile(pwd,'mpcCustomSolver.m'); copyfile(src,dest,'f')
Симулируйте управление с обратной связью двигателя постоянного сервопривода и сохраните вход и выход объекта.
sim(mdl) uDantzigSim = u; yDantzigSim = y;
-->Converting model to discrete time. Assuming no disturbance added to measured output channel #1. -->The "Model.Noise" property of the "mpc" object is empty. Assuming white noise on each measured output channel.
Чтобы запустить оставшийся пример, требуется продукт Simulink Coder.
if ~mpcchecktoolboxinstalled('simulinkcoder') disp('Simulink(R) Coder(TM) is required to run this example.') return end
Чтобы сгенерировать код из блока MPC Controller, который использует пользовательский решатель QP, включите пользовательский решатель для генерации кода, опции в контроллере MPC.
mpcobj.Optimizer.CustomSolverCodeGen = true;
Вы также должны предоставить функцию MATLAB ®, которая удовлетворяет всем следующим требованиям:
Имя функции должно быть mpcCustomSolverCodeGen.
Входные и выходные аргументы должны совпадать с аргументами в файле шаблона.
Функция должна находиться в пути MATLAB.
В этом примере используйте тот же пользовательский решатель, определенный в mpcCustomSolverCodeGen_TemplateEML.txt. Сохраните функцию в рабочей папке следующим mpcCustomSolverCodeGen.m.
src = which('mpcCustomSolverCodeGen_TemplateEML.txt'); dest = fullfile(pwd,'mpcCustomSolverCodeGen.m'); copyfile(src,dest,'f')
Проверьте сохраненные mpcCustomSolverCodeGen.m файл.
function [x, status] = mpcCustomSolverCodeGen(H, f, A, b, x0) %#codegen % mpcCustomSolverCodeGen allows the user to specify a custom (QP) solver % written in MATLAB to be used by MPC controller in code generation. % % Workflow: % (1) Copy this template file to your work folder and rename it to % "mpcCustomSolverCodeGen.m". The work folder must be on the path. % (2) Modify the "mpcCustomSolverCodeGen.m" to use your solver. % Note that your solver must use only fixed-size data. % (3) Set "mpcobj.Optimizer.CustomSolverCodeGen = true" to tell the MPC % controller to use the solver in code generation. % To generate code: % In MATLAB, use "codegen" command with "mpcmoveCodeGeneration" (require MATLAB Coder) % In Simulink, generate code with MPC and Adaptive MPC blocks % % To use this solver for simulation in MATLAB and Simulink, you need to: % (1) Copy "mpcCustomSolver.txt" template file to your work folder and % rename it to "mpcCustomSolver.m". The work folder must be on the path. % (2) Modify the "mpcCustomSolver.m" to use your solver. % (3) Set "mpcobj.Optimizer.CustomSolver = true" to tell the MPC % controller to use the solver in simulation. % % The MPC QP problem is defined as follows: % % min J(x) = 0.5*x'*H*x + f'*x, s.t. A*x >= b. % % Inputs (provided by MPC controller at run-time): % H: a n-by-n Hessian matrix, which is symmetric and positive definite. % f: a n-by-1 column vector. % A: a m-by-n matrix of inequality constraint coefficients. % b: a m-by-1 vector of the right-hand side of inequality constraints. % x0: a n-by-1 vector of the initial guess of the optimal solution. % % Outputs (sent back to MPC controller at run-time): % x: must be a n-by-1 vector of optimal solution. % status: must be an integer of: % positive value: number of iterations used in computation % 0: maximum number of iterations reached % -1: QP is infeasible % -2: Failed to find a solution due to other reasons % Note that: % (1) When solver fails to find an optimal solution (status<=0), "x" % still needs to be returned. % (2) To use sub-optimal QP solution in MPC, return the sub-optimal "x" % with "status = 0". In addition, you need to set % "mpcobj.Optimizer.UseSuboptimalSolution = true" in MPC controller. % % DO NOT CHANGE LINES ABOVE % This template implements a showcase QP solver using "Dantzig" algorithm % by G. B. Dantzig, A. Orden, and P. Wolfe, "The generalized simplex method % for minimizing a linear form under linear inequality constraints", % Pacific J. of Mathematics, 5:183–195, 1955. % % User is expected to modify this template and plug in own custom QP solver % that replaces the "Dantzig" algorithm. ZERO = zeros('like',H); ONE = ones('like',H); % xmin is a constant term that adds to the initial basis because "dantzig" % requires positive optimization variables. A fixed "xmin" does not work % for all MPC problems. xmin = -1e3*ones(size(f(:)))*ONE; maxiter = 200*ONE; nvar = length(f); ncon = length(b); a = -H*xmin(:); H = H\eye(nvar); rhsc = A*xmin(:) - b(:); rhsa = a-f(:); TAB = -[H H*A';A*H A*H*A']; basisi = [H*rhsa; rhsc + A*H*rhsa]; ibi = -(1:nvar+ncon)'*ONE; ili = -ibi*ONE; %% Call EML function "qpdantzg" [basis,ib,il,iter] = qpdantzg(TAB,basisi,ibi,ili,maxiter); %#ok<ASGLU> %% status if iter > maxiter status = ZERO; elseif iter < ZERO status = -ONE; else status = iter; end %% optimal variable x = zeros(nvar,1,'like',H); for j = 1:nvar if il(j) <= ZERO x(j) = xmin(j); else x(j) = basis(il(j))+xmin(j); end end
Сгенерируйте исполняемый код из модели Simulink с помощью slbuild команда от Simulink Coder.
slbuild(mdl)
### Starting build procedure for: mpc_customQPcodegen -->Converting model to discrete time. Assuming no disturbance added to measured output channel #1. -->The "Model.Noise" property of the "mpc" object is empty. Assuming white noise on each measured output channel. ### Successful completion of build procedure for: mpc_customQPcodegen Build Summary Top model targets built: Model Action Rebuild Reason ==================================================================================================== mpc_customQPcodegen Code generated and compiled Code generation information file does not exist. 1 of 1 models built (0 models already up to date) Build duration: 0h 0m 22.665s
В системе Windows после концов процесса сборки программное обеспечение добавляет исполняемый файл mpc_customQPcodegen.exe в рабочую папку.
Запустите исполняемый файл. После того, как исполняемый файл успешно завершится (status = 0), программное обеспечение добавляет файл данных mpc_customQPcodegen.mat в рабочую папку. Загрузите файл данных в рабочее пространство MATLAB и получите входные и выходные сигналы объекта, сгенерированные исполняемым файлом.
if ispc status = system(mdl); load(mdl) uDantzigCodeGen = u; yDantzigCodeGen = y; else disp('The example only runs the executable on Windows system.'); end
The example only runs the executable on Windows system.
Сравните входные и выходные сигналы объекта из всех симуляций.
if ispc figure subplot(2,1,1) plot(u.time,uKWIK.signals.values,u.time,uDantzigSim.signals.values,... '+',u.time,uDantzigCodeGen.signals.values,'o') subplot(2,1,2) plot(y.time,yKWIK.signals.values,y.time,yDantzigSim.signals.values,... '+',y.time,yDantzigCodeGen.signals.values,'o') legend('KWIK','Dantzig Simu','Dantzig CodeGen','Location','northwest') else figure subplot(2,1,1) plot(u.time,uKWIK.signals.values,u.time,uDantzigSim.signals.values,'+') subplot(2,1,2) plot(y.time,yKWIK.signals.values,y.time,yDantzigSim.signals.values,'+') legend('KWIK','Dantzig Simu','Location','northwest') end
Сигналы от всех симуляций идентичны.
[1] Bemporad, A. and Mosca, E. «Выполнение жестких ограничений в неопределенных линейных системах путем ссылки». Automatica, Vol. 34, Number 4, pp. 451-461, 1998.
bdclose(mdl)