Приложение MPC с онлайн-мониторингом

Установка, используемая в этом примере, представляет собой систему с одним входом и одним выходом с жесткими ограничениями как на управляемую переменную (MV), так и на управляемый выход (OV). Целью управления является удержание OV на уставке 0. К OV добавляется неизмеренное нарушение нагрузки. Это возмущение изначально представляет собой увеличение клина. Отклик контроллера в конечном итоге насыщает СН на его жестком пределе. Как только происходит насыщение, контроллер больше ничего не может сделать, и нарушение в конечном итоге приводит OV выше заданного жесткого верхнего предела. Когда контроллер прогнозирует, что невозможно заставить OV ниже этого верхнего предела, QP во время выполнения становится неосуществимым.

Определите установку как систему SISO первого порядка с единичным коэффициентом усиления.

Plant = tf(1,[2 1]);

Определите неизмеренное нарушение нагрузки. Сигнал нарастает от 0 кому 2 между 1 и 3 секунды, затем наклоняется обратно вниз от 2 кому 0 между 3 и 5 секунд.

LoadDist = [0 0; 1 0; 3 2; 5 0; 7 0];

Проектирование контроллера MPC

Создайте объект MPC с использованием модели тестового завода. Выбранный интервал регулирования составляет примерно одну десятую доминирующей постоянной времени растения.

Ts = 0.2;
Obj = mpc(Plant, Ts);

-->The "PredictionHorizon" property of "mpc" object is empty. Trying PredictionHorizon = 10.
-->The "ControlHorizon" property of the "mpc" object is empty. Assuming 2.
-->The "Weights.ManipulatedVariables" property of "mpc" object is empty. Assuming default 0.00000.
-->The "Weights.ManipulatedVariablesRate" property of "mpc" object is empty. Assuming default 0.10000.
-->The "Weights.OutputVariables" property of "mpc" object is empty. Assuming default 1.00000.

Определите жесткие ограничения на вход (СН) и выход (OV) установки. По умолчанию все ограничения MV являются жесткими, а ограничения OV - мягкими.

Obj.MV.Min = -0.5;
Obj.MV.Max =  1;
Obj.OV.Min = -1;
Obj.OV.Max =  1;
Obj.OV.MinECR = 0; % change OV lower limit from soft to hard
Obj.OV.MaxECR = 0; % change OV upper limit from soft to hard

Как правило, жесткие ограничения OV не поощряются и используются здесь только для иллюстрации того, как обнаружить неосуществимое QP. Жесткие ограничения OV делают несходимость вероятной, и в этом случае необходим план управления резервным копированием. Этот пример не включает план резервного копирования. Однако, как показано в моделировании, действие по умолчанию замораживания одиночного СН является наилучшим ответом в этом простом случае.

Моделирование с использованием Simulink со встроенным решателем QP

Для выполнения этого примера необходимы Simulink и панель инструментов оптимизации.

if ~mpcchecktoolboxinstalled('simulink')
    disp('Simulink is required to run this example.')
    return
end
if ~mpcchecktoolboxinstalled('optim')
    disp('The Optimization Toolbox is required to run this example.')
    return
end

Создайте систему управления в модели Simulink и включите qp.status outport путем выбора параметра Optimization status блока MPC Controller. Просмотр времени выполнения qp.status значение в области состояния контроллера.

mdl = 'mpc_onlinemonitoring';
open_system(mdl)

Смоделировать отклик с замкнутым контуром с помощью решателя QP панели управления прогнозом модели по умолчанию.

open_system([mdl '/Controller Status'])
open_system([mdl '/Response'])
sim(mdl)

-->Converting the "Model.Plant" property of "mpc" object to state-space.
-->Converting model to discrete time.
-->Assuming output disturbance added to measured output channel #1 is integrated white noise.
-->The "Model.Noise" property of the "mpc" object is empty. Assuming white noise on each measured output channel.

Объяснение ответа с замкнутым контуром

Как показано в объеме ответа, через 1,4 секунды возрастающее возмущение заставляет СН насыщаться на нижней границе -0,5, что является решением QP в этих условиях (потому что контроллер пытается удерживать OV на своей уставке 0).

OV продолжает увеличиваться из-за нарушения клина и через 2,2 секунды превышает указанную жесткую верхнюю границу 1.0. Поскольку QP формулируется в терминах прогнозируемых выходных данных, контроллер все еще прогнозирует, что он может вернуть OV ниже 1,0 в следующем движении, и поэтому проблема QP все еще выполнима.

Наконец, при t = 3,2 секунды контроллер прогнозирует, что он больше не может перемещать OV ниже 1,0 в течение следующего интервала управления, и проблема QP становится неосуществимой и qp.status в настоящее время изменяется на -1.

Через три секунды возмущение уменьшается. Через 3,8 секунды QP снова становится возможным. Однако OV все еще значительно выше своей уставки, и MV остается насыщенным до 5,4 секунды, когда решение QP должно увеличить MV, как показано. С этого момента СН не насыщается, и контроллер может вернуть OV обратно к своей уставке.

Когда QP осуществим, встроенный решатель находит решение в трех итерациях или менее.

Моделирование с помощью пользовательского решателя QP

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

Obj.Optimizer.CustomSolver = true;

Необходимо также предоставить функцию MATLAB ®, которая удовлетворяет всем следующим требованиям:

Для этого примера используйте пользовательский решатель, определенный в mpcCustomSolver.txt, которая использует quadprog на панели инструментов оптимизации в качестве пользовательского решателя QP. Чтобы внедрить собственный решатель QP, измените этот файл.

Сохраните функцию в рабочей папке как .m файл.

src = which('mpcCustomSolver.txt');
dest = fullfile(pwd,'mpcCustomSolver.m');
copyfile(src,dest,'f');

Просмотреть сохраненные mpcCustomSolver.m файл.

function [x, status] = mpcCustomSolver(H, f, A, b, x0)
% mpcCustomSolver allows user to specify a custom quadratic programming
% (QP) solver to solve the QP problem formulated by MPC controller.  When
% the "mpcobj.Optimizer.CustomSolver" property is set true, instead of
% using the built-in QP solver, MPC controller will now use the customer QP
% solver defined in this function for simulations in MATLAB and Simulink.
%
% The MPC QP problem is defined as follows:
%   Find an optimal solution, x, that minimizes the quadratic objective
%   function, J = 0.5*x'*H*x + f'*x, subject to linear inequality
%   constraints, 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 (fed back to MPC controller at run-time):
%       x: must be a n-by-1 vector of optimal solution. 
%  status: must be an finite 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 even if solver failed to find an optimal solution, "x" must be
% returned as a n-by-1 vector (i.e. set it to the initial guess x0)
%
% DO NOT CHANGE LINES ABOVE

% The following code is an example of how to implement the custom QP solver
% in this function.  It requires Optimization Toolbox to run.

% Define QUADPROG options and turn off display of optimization results in
% Command window.
options = optimoptions('quadprog');
options.Display = 'none';  
% By definition, constraints required by "quadprog" solver is defined as
% A*x <= b.  However, in our MPC QP problem, the constraints are defined as
% A*x >= b.  Therefore, we need to implement some conversion here:
A_custom = -A;
b_custom = -b;
% Compute the QP's optimal solution.  Note that the default algorithm used
% by "quadprog" ('interior-point-convex') ignores x0.  "x0" is used here as
% an input argument for illustration only.
H = (H+H')/2; % ensure Hessian is symmetric
[x, ~, Flag, Output] = quadprog(H, f, A_custom, b_custom, [], [], [], [], x0, options);
% Converts the "flag" output to "status" required by the MPC controller.
switch Flag
    case 1
        status = Output.iterations;
    case 0
        status = 0;
    case -2
        status = -1;
    otherwise
        status = -2;
end
% Always return a non-empty x of the correct size.  When the solver fails,
% one convenient solution is to set x to the initial guess.
if status <= 0
    x = x0;
end

Повторите моделирование.

set_param([mdl '/Controller Status'],'ymax','10');
sim(mdl)

-->Converting the "Model.Plant" property of "mpc" object to state-space.
-->Converting model to discrete time.
-->Assuming output disturbance added to measured output channel #1 is integrated white noise.
-->The "Model.Noise" property of the "mpc" object is empty. Assuming white noise on each measured output channel.

Входные и выходные сигналы установки идентичны сигналам, полученным с помощью встроенного решателя панели управления прогнозом модели, но qp.status показывает, что quadprog не занимает одинаковое количество итераций, чтобы найти решение. Однако он обнаруживает тот же период времени несходимости.

bdclose(mdl);