Модель сервопривода DC

Модель dc-сервопривода является линейной динамической системой, описанной в [1]. plant модель в пространстве состояний непрерывного времени двигателя. tau максимальный допустимый крутящий момент, который вы используете в качестве выходного ограничения.

[plant,tau] = mpcmotormodel;

Спроектируйте контроллер MPC

Объект имеет вход того, моторное входное напряжение. Диспетчер 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 со встроенным решателем QP

Чтобы запустить остающийся пример, 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;

Симулируйте в Simulink с пользовательским решателем QP

Чтобы исследовать, как пользовательский решатель ведет себя при тех же условиях, включите пользовательский решатель в контроллере MPC.

mpcobj.Optimizer.CustomSolver = true;

Необходимо также предусмотреть функцию 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.

Сгенерируйте код с пользовательским решателем QP

Чтобы запустить остающийся пример, продукт 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_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 Embedded 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 Embedded MATLAB 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 (require Simuink Coder products)
%
% 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 с помощью rtwbuild команда от Simulink Coder.

rtwbuild(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

В системе 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)