Fitting operating point for reduced order model VERY inaccurate

[yonghoonlee]

@danielrherber and @AthulKrishnaSundarrajan

I tried several different ways to fit the state variables at the operating point, including least square fitting without scaling, with scaling, and mldivide calculations. None of them gives good solution.

Maybe (as Dan previously mentioned) fit the state variable to a few important output variables could be more promising, especially because of the fact that output variables have varied range of scales.

Using the full 102-state model, (servo output power) \approx (torque)*(angspeed) \approx 15MW. The errors were at most a few hundred kW. However, the reduced model with fitted state values at the operating point exhibits about errors around 2MW. This level of accuracy could make this model not very useful.

Let's discuss this more in depth once @AthulKrishnaSundarrajan setup his development environment. Please let me know. Thank you.

[danielrherber]

Have you tried using the Kalman decomposition of (A,B,C)?

[sysr,u] = minreal(sys,tol) returns, for state-space model sys, an orthogonal matrix U such that (U*A*U',U*B,C*U') is a Kalman decomposition of (A,B,C)

https://www.mathworks.com/help/control/ref/minreal.html

[danielrherber]

Ignore my previous comment above. I believe I have found the mathematically rigorous way to do this using balreal. I have two tests below. Essentially, we need to use the transformation matrix on the initial states, then add the contribution of eliminated states to the new initial condition. The assumption is that those states don't change the output value much (hence they were eliminated).

[danielrherber]

close all; clear; clc

testnum = 2;

switch testnum
    %----------------------------------------------------------------------
    case 1
    s = tf('s');
    G11 = 1/(s+1);
    G12 = 2/(s+1);
    G21 = -1/(s+1);
    G22 = -1/(s+1);
    R1 = 5*(s+1)/s;
    tf = -G12*G21*R1 /(1+R1*G11);      
    sys = ss(tf);
    %----------------------------------------------------------------------
    case 2
    rng(3635)
    ny = 100; nu = 1; no = 2;
    A = full(sprand(ny,ny,0.1));
    B = full(sprand(ny,nu,0.1));
    C = full(sprand(no,ny,0.1));
    D = [];
    sys = ss(A,B,C,D);
        
    %----------------------------------------------------------------------
end

[sysbal,g,T,Te] = balreal(sys);  % Compute balanced realization
elim = (g<1e-8);         % Small entries of g are negligible states
sysred = modred(sysbal,elim); % Remove negligible states

% initial conditions
X0 = 10.^(3*rand(length(g),1));

% original output
disp("original initial output value")
disp(sys.C*X0)

% balanced output
disp("balanced initial output value")
disp(sysbal.C*T*X0)

% extract parts of T (eliminated and kept states)
Tk = T; Tk(elim,:) = []; Tk(:,elim) = [];
Te = T; Te(~elim,:) = []; Te(:,~elim) = [];

% compute offset from eliminated states initial values
CT = sysbal.C*T;
Ce = CT(:,elim);
X0e = X0(elim);


% reduced output
disp("reduced initial output value")
disp(CT(:,~elim)*X0(~elim) + Ce*X0e) % <- key equation, note the offset

%--- try more initial states
disp(" ")

% number of tests
ntest = 10;

% go through each test
for k = 1:ntest
    
   % new random initial states
   X0 = 10.^(3*rand(length(g),1));
   
   % keep eliminated states fixed
   X0(elim) = X0e;
   
   % original initial output value 
   y0 = sys.C*X0;
   
   % reduced initial output value
   yred = CT(:,~elim)*X0(~elim) + Ce*X0e;
   % yred = sysred.C*(Tk*X0(~elim)) + Ce*X0e;
   
   % display the relative error in %
   disp(100*[(y0-yred)./(y0 + 1e-14)]')
   
end

[yonghoonlee]

@danielrherber I did not check how much more accurate, but Kalman decomposition is much much more better than all the other methods. I just saw your last comment after I implemented Kalman decomposition, so I will try this as well.

[yonghoonlee]

However, it seems that if we use minreal, the minreal internally uses Kalman decomposition to get the reduced order model. So, Kalman decomposition-based method should give us quite accurate value, in case reduced order model is not wrong or inaccurate. I will update once I get more test results next week.

[danielrherber]

@yonghoonlee Great! The relative error in the tests was less than 1% using many random states which I believe implies that the output error in the reduced model would never be greater than 1% (or so, numbers approximate for this test).

[danielrherber]

Sorry, one more update. I fixed something in the key equation. This approach is now is O(1e-10) accuracy. The assumption that the eliminated states don't change is key and that the current states are accurate (with respect to the original model).

[danielrherber]

@yonghoonlee I am not sure if you followed up with this approach but tomorrow we will present some validation of the reduced linear model using balreal and the difference in the outputs between the models.

[danielrherber]

Work-in-progress script for validating the linear model and the reduced linear model.
TEST_Balanced_Realization.zip