Hi,
I am currently working on a control system for an over-actuated aircraft with adaptive on-board actuator dynamics model. This actuator dynamics model is used within the control allocation optimizer. To estimate the actuator dynamics online, I use a recursive least squares (RLS) algorithm with adaptive forgetting factor, based on the information that is available in the input signal (to omit persistence of excitation requirements).
My control allocation optimizer can only take first order actuator dynamics, because my control effectiveness continuity is limited and I want to omit the need for higher order output within the control scheme. The RLS estimator however, estimates the parameters of a second order actuator dynamics model. To use this information within my control allocator, I need to find the best fitting first order lag (gain+time constant) for the second order dynamics estimate (estimating the parameters of a first order lag directly won't work, as it gets really confused about overshoot and stuff).
Currently, my actuator dynamics parameter update model code looks like this, with pk the parameter estimates and Pk the covariance at timestep k. In the second for loop, I translate the parameter estimates of the second order difference equation into a discrete state space model, after which I try to reduce the model to a first order lag. For this I use Task-Based Model Order Reduction Workflow - MATLAB & Simulink - MathWorks Benelux. However, because this code is within my Simulink Matlab fcn block, I need to call coder.extrinsic, which makes my simulation rather slow. Is there anybody who can help me to make the simulation go a bit faster again. Either by omitting the task-based model reduction for some other method, or any other way?
Regards,
Noah
function [Pk, pk] = fcn(P0,hk,sig,yk,lam_min,d_on,p0)
% Function to estimate the parameters a, b, c, d of
% d_k+1 = a*d_k + b*d_k-1 + c*u_k + d*u_k-1
% Define persistent variables
persistent p P
if isempty(p)
p=p0;
end
if isempty(P)
P=P0;
end
% dimensions hk [4x13], yk [13x4], pk [4x13], Pk [4x13x4]
for i = d_on' % for each actuator
Pi = squeeze(P(:,i,:)); % [4x4]
hi = hk(i,:)'; % [4x1]
yi = yk(i); % [1]
pi = p(:,i); % [4x1]
% Compute least squares gain L and adaptive forgetting factor lambda
err = yi-hi'*pi; % [1]
L = Pi*hi*inv(hi'*Pi*hi+1); % [4x1]
lam = 1-(1-hi'*L)*err^2/sig;
if lam < lam_min
lam = lam_min;
end
% Determine new parameter estimates for i-th actuator
p(:,i) = pi+L*err; % [4x1]
% Determine new covariance matrix for i-th actuator
P(:,i,:) = (eye(4)-L*hi')*Pi/lam; % [4x4]
end
% Put into control canonical state space form and reduce to first order lag
pk=zeros(2,13);
coder.extrinsic('ss', 'reducespec', 'getrom', 'ssdata', 'process')
for i = d_on'
% Extract parameter estimates
a = p(1,i);
b = p(2,i);
c = p(3,i);
d = p(4,i);
% Define state space model (control canonicle form)
PHI = [a, b; 1, 0];
GAM = [1; 0];
PSI = [c, d];
sys = ss(PHI,GAM,PSI,0,0.01);
% Simplify to first order actuator dynamics model estimates
R = reducespec(sys,"balanced");
R = process(R);
rsys = getrom(R,Order=1);
[p1,p2,p3,~] = ssdata(rsys);
pk(1,i) = p1;
pk(2,i) = p3*p2;
end
Pk = P;
