How to use "StabilityThreshold" estimation option in greyest function?

Question

준호 2024-10-21，7:55

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2161030-how-to-use-stabilitythreshold-estimation-option-in-greyest-function

回答： Divyanshu 2024-10-28，4:01

Hello,

To predict the lateral equations of the H4 Hercules using a greybox model, I defined the model function as follows and set the initial values for the parameters before executing it. (Since I am estimating a linear model, I extracted the data from the trim state, so the initial values for all state variables and control inputs are set to 0.)

function [A,B,C,D] = h4_lat_dynamics(theta, ~)
    a11 = theta(1); a12 = theta(2); a13 = theta(3); a14 = theta(4);
    a21 = theta(5); a22 = theta(6); a23 = theta(7);
    a31 = theta(8); a32 = theta(9); a33 = theta(10); a36 = theta(11); 
    a37 = theta(12); a38 = theta(13); a39 = theta(14); a310 = theta(15); 
    a311 = theta(16); a312 = theta(17); a313 = theta(18); a66 = theta(19);
    a77 = theta(20); a88 = theta(21); a99 = theta(22); a1010 = theta(23); 
    a1111 = theta(24); a1212 = theta(25); a1313 = theta(26);
    b12 = theta(27); b21 = theta(28); b22 = theta(29); b31 = theta(30); 
    b32 = theta(31); b66 = theta(32); b77 = theta(33); b88 = theta(34); 
    b99 = theta(35); b1010 = theta(36); b1111 = theta(37); 
    b1212 = theta(38); b1313 = theta(39);
    A = [a11 a12 a13 a14 0 0 0 0 0 0 0 0 0;
         a21 a22 a23  0  0 0 0 0 0 0 0 0 0;
         a31 a32 a33  0  0 a36 a37 a38 a39 a310 a311 a312 a313;
          0   1   0   0  0 0 0 0 0 0 0 0 0;
          0   0   1   0  0 0 0 0 0 0 0 0 0;
          0   0   0   0  0 a66 0 0 0 0 0 0 0;
          0   0   0   0  0 0 a77 0 0 0 0 0 0;
          0   0   0   0  0 0 0 a88 0 0 0 0 0;
          0   0   0   0  0 0 0 0 a99 0 0 0 0;
          0   0   0   0  0 0 0 0 0 a1010 0 0 0;
          0   0   0   0  0 0 0 0 0 0 a1111 0 0;
          0   0   0   0  0 0 0 0 0 0 0 a1212 0;
          0   0   0   0  0 0 0 0 0 0 0 0 a1313];
    
    B = [0   b12 0 0 0 0 0 0 0 0;
         b21 b22 0 0 0 0 0 0 0 0;
         b31 b32 0 0 0 0 0 0 0 0;
         0   0   0 0 0 0 0 0 0 0;
         0   0   0 0 0 0 0 0 0 0;
         0   0   b66 0 0 0 0 0 0 0;
         0   0   0 b77 0 0 0 0 0 0;
         0   0   0 0 b88 0 0 0 0 0;
         0   0   0 0 0 b99 0 0 0 0;
         0   0   0 0 0 0 b1010 0 0 0;
         0   0   0 0 0 0 0 b1111 0 0;
         0   0   0 0 0 0 0 0 b1212 0;
         0   0   0 0 0 0 0 0 0 b1313];
    C = eye(size(A,1));
    D = zeros(size(C,1),size(B,2));
    
    sys = c2d(ss(A,B,C,D),1/10);
    
    A = sys.A;
    B = sys.B;
    C = sys.C;
    D = sys.D;
end

z = iddata(states_delta_lat,controls_delta_lat,0.1,'Name','Data');
z.InputName = {'Aileron', 'Rudder', 'First-engine throttle', 'Second-engine throttle', 'Third-engine throttle', 'Fourth-engine throttle',...
'Fifth-engine throttle', 'Sixth-engine throttle', 'Seventh-engine throttle', 'Eighth-engine throttle'};
z.InputUnit = {'normalized', 'normalized', 'normalized', 'normalized', 'normalized', 'normalized', 'normalized', 'normalized',...
    'normalized', 'normalized'};
z.OutputName = {'lateral velocity', 'roll-rate', 'yaw-rate', 'roll angle', 'yaw angle', 'First-engine thrust', 'Second-engine thrust',...
    'Thrid-engine thrust', 'Fourth-engine thrust', 'Fifth-engine thrust', 'Sixth-engine thrust', 'Seventh-engine thrust', 'Eighth-engine thrust'};
z.OutputUnit = {'ft/s', 'deg/s', 'deg/s', 'deg', 'deg', 'lbf', 'lbf', 'lbf', 'lbf','lbf', 'lbf', 'lbf', 'lbf'};
z.Tstart = 0;
z.TimeUnit = 's';
init = {[-0.005; 0.002; -0.998; 0.0829; -1.5800; -0.5600; 0.2560; 0.3588; -0.0980; -0.1600; -0.008; -0.006; -0.004; -0.002; 0.002; 0.004; 0.006; 0.008;...
    1.000; 1.000; 1.000; 1.000; 1.000; 1.000; 1.000; 1.000; 0.0250; 0.1500; 0.1900; 0.0260; -0.3660; 25.00; 25.00; 25.00; 25.00; 25.00; 25.00; 25.00; 25.00]};
aux = {};
init_sys = idgrey('h4_lat_dynamics', init, 'c', aux, 0);
opt = greyestOptions;
opt.DisturbanceModel = 'none';
opt.EnforceStability = true;
h4_lat_est = greyest(z,init_sys);
compare(z,h4_lat_est,'-b');

However, after running the code, I received the following warning, and upon comparing the original data with the estimated data, I found that the model was highly unstable.

Warning: The initial model has an unstable predictor (some eigenvalues of (A-K*C) are outside the stability threshold). You can use the "init" command to create a stable predictor. Use the "StabilityThreshold" estimation option to set custom thresholds.