В этом примере показано, как симулировать и сгенерировать код для прогнозирующего контроллера модели, который использует пользовательский решатель квадратичного программирования (QP). Объект для этого примера является dc-сервоприводом в Simulink®.
Модель dc-сервопривода является линейной динамической системой, описанной в [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, которая симулирует управление с обратной связью dc-сервопривода с помощью контроллера 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')
Симулируйте управление с обратной связью dc-сервопривода и сохраните ввод и вывод объекта.
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 28.016s
В системе 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. и Моска, E. "Выполняя трудные ограничения в неопределенных линейных системах ссылочным управлением". Automatica, Издание 34, Номер 4, стр 451-461, 1998.
bdclose(mdl)