Идентификация модели ARX с ограничениями линейного равенства и стабильностью

Question

Из входных и выходных данных я должен идентифицировать с наименьшими квадратами модель ARX

y(n) = -a1 y(n-1) -....- aN y(nN) + b1 x(n-1)+...+bM x(нМ)

который должен быть устойчивым (все его полюсы находятся внутри единичного круга) и в то же время параметры

a1,...,aN, b1,...,bM

должны удовлетворять некоторым ограничениям линейного равенства. Например, мне нужно статическое усиление, равное 1, это означает,

-a1-a2...-aN+b1+b2...+bM =1

В матричной форме я могу написать ограничение как:

Тета = [a1 a2... aN b1 b2... bM] ';

Aeq = [-он (1,N) единиц (1, M)];

beq = 1;

Aeq * тета = BEQ;

С помощью функции Matlab ниже arx_stable Я могу идентифицировать конюшню ARX используя функцию Matlab 'ARX' ​​для установки параметра 'Focus' ='stability' (но без равенства линейных ограничений)

С функцией ниже arx_constr Я могу определить ARX с линейным ограничением равенства Aeq beq (но без устойчивости)

Как я могу получить стабильность и уважение ограничений??? Спасибо,

Симона

function theta = arx_stable(N,M,t,input_data,output_data)
% input_data and output_data are column vectors
% t is the time vector
% data structure generation:
Ts = t(2)-t(1); % it is the sampling time
data = iddata(output_data,input_data,Ts);
n_k=0; % it is the Input-output delay
ARX_model = arx(data,[N M n_k]  , 'Focus','stability');  
% 'Focus','stability' in this way the identified ARX is stable
a=ARX_model.a;
a=a(2:end); % I remove the first a0=1
b=ARX_model.b;
theta=[a b];
end

function theta = arx_constr(N,M,t,input_data,output_data,Aeq,beq)
% input_data and ouput_data are row vectors
% t is the time vector
% X and y construction:
n_order=max([N M]);
X=[];  % X is the regression matrix
for nn = n_order+1:length(input_data)
    X(nn,:)=[-output_data((nn-1):-1:(nn-N)) input_data((nn-1):-1:(nn-M))];
end
X=X((n_order+1):end,:); % it cancels the first n_order rows that are full of zeros
y=output_data(n_order+1:end)' ;   %  y is a column vector
options_lsqlin=optimset('Algorithm','active-set','LargeScale','off');   
% this set is necessary to use equality constraints
theta=lsqlin(X,y,[],[],Aeq,beq,[],[],[],options_lsqlin);
end