Well-identified fitted process model does not behave like data on simulink

Hello, and thanks for taking your time to read my question.

I have a set of experiments that I am trying to identify using the system identification tool, and I succesfully obtained a process model (2nd order transfer function, 2 inputs 1 output + delay each) and a state-space model (order 2). The fits in the identification tool (model output window) are almost perfect but when I try to simulate the same data with the process model (idmodel) block in simulink the output is completely messed up. I used the exact same simulink model with the state-space fitted model and the output is congruent with the toolbox model output window. I think the reason is because I can't set initial values for transfer function identified models like I can do for state-space models .

Is there a way to set initial values for transfer function models on simulink? I already tried converting the process model to a state space model with ss() and use a state space block but it didn't work, I assume because in the model options I selected a first order disturbance term. Also I tried fitting the model with zero initial conditions but the results were the same (bad). I posted this question here as well.

My objective is to simulate the identified process model just as it appears in the system identification tool (model output) window, maybe the current approach is not the correct one so it would be absolute beneficial to me if anyone knows how to do so and is kind enough to set me on the right path.

Figure 1. Model output window. Yellow is SS model, blue is process model

Figure 2. Process model parameter window. Fit options are as displayed, everything else is on default values.

Figure 3. Simulation of process model with simulink. It can be observed the fit (blue) does not correspond to data (yellow)

Figure 4. Simulation of state space model with simulink.It is coherent with the model output window of figure 1

Figure 5. Simulink idmodel block parameters. Left: state space model, it allows to enter initial conditions. Right: process model. Can't set initial conditions.

[julian navas]

Update

I can compare the models and the data using [y,fit, x0]=compare(sys, data) but still can't figure out how to do this in simulink. Please any help is greatly welcome.

Answer 1

[Rajiv Singh]

You are almost there. Convert the model into state-space form and use it for simulation. For initial conditions, you will need to run FINDSTATES command on the dataset whose output you are trying to match. See:

https://www.mathworks.com/help/ident/ug/estimate-initial-conditions.html

[julian navas]

Hello and thanks for reviewing my question. Maybe I didn't put it into words correctly but I already did what you suggest, please read the 2nd paragraph and the text in the figure 5 initial states parameter. When I convert the model to ss I use ss(sys,'augmented') since it has the noise component of the previous model as a new state (as far as I understand) resulting in a 5 state matrix form. I then use findstates. It does not give good results. The approach you suggest involves converting the system using idss, I didn't know of this method, will give it a try and see if (hope) it works. Thank you.

[julian navas]

I tried your answer and got mixed results. Let me explain: I run this to convert my idproc model to ss and to estimate initial states.

model=idss(T3s_2d);
findstates(model,run_data_s{8})
ans =

         0
         0
  186.2812
  -34.2960
  106.0995
  
findstates(T3s_2d,run_data_s{8})

ans =

         0
         0
  186.2812
  -34.2960
  106.0995
  

Although the SS matrices are not in the canonical form of the other models, it doesn't matter, the initial conditions are the same. I use compare to see what their output is like on the model output window of the system identification toolbox

  compare(run_data_s{8}, T3s_2d, model)
  

Excellent. Both appear with the same name because that's the internal name in the model properties even after conversion. Now I try it in simulink but it doesn't work as expected.

I thought that maybe I was missing noise information, then checked the 'Add noise' button in the idmodel parameters but it got much much more worse.

This is better than what I initially had, so it's progress, but still not quite there yet. Do you have more ideas on how to handle this strange (for me) behaviour?

[julian navas]

Update

Since the model output window shows good results I checked the initial conditions x0 returned by compare:

[y,fit,x0]=compare(run_data_s{8}, model, T3s_2d);
cell2mat(x0)
ans =

   1.0e+04 *

    3.6829
   -1.9391
    0.0379
    0.1520
    0.0569
    3.6829
   -1.9391
    0.0379
    0.1520
    0.0569
    

Both, again, have the same set of initial conditions, however this x0 set is very different to the one found with findstates. I used this set of initial conditions on the simulink block and this time it worked

[y,fit,x0]=compare(run_data_s{8}, model);

I finally accomplished my goal but still have questions. Why does findstates() and compare() give different set of initial conditions?

In any case, thank you so mucho for pointing me in the right direction. You're awesome.

[Rajiv Singh]

can you print the model equation? I want to see number of poles/zeros and delays for each I/O pair

[julian navas]

Of course, here it is

>> model=idss(T3s_2d)

model =
  Continuous-time identified state-space model:
      dx/dt = A x(t) + B u(t) + K e(t)
       y(t) = C x(t) + D u(t) + e(t)
 
  A = 
               x1          x2          x3          x4          x5
   x1   -0.001739  -0.0003834           0           0           0
   x2   0.0002441           0           0           0           0
   x3           0           0    -0.02643  -0.0006169           0
   x4           0           0   0.0009766           0           0
   x5           0           0           0           0   -0.002102
 
  B = 
        T_{pv}   P_{pv}
   x1  0.01563        0
   x2        0        0
   x3        0     0.25
   x4        0        0
   x5        0        0
 
  C = 
             x1       x2       x3       x4       x5
   LOD        0  0.01396        0   0.1303   0.1476
 
  D = 
        T_{pv}  P_{pv}
   LOD       0       0
 
  K = 
        LOD
   x1     0
   x2     0
   x3     0
   x4     0
   x5  0.25
 
  Input delays (seconds): 630  8.79 
 
Name: T3s_2d
Parameterization:
   FREE form (all coefficients in A, B, C free).
   Feedthrough: none
   Disturbance component: estimate
   Number of free coefficients: 45
   Use "idssdata", "getpvec", "getcov" for parameters and their uncertainties.

Status:                                 
Created by conversion from idproc model.
>> T3s_2d

T3s_2d =
Process model with 2 inputs: y = G11(s)u1 + G12(s)u2         
  From input "T_{pv}" to output "LOD":                       
                  Kp                                         
  G11(s) = ----------------- * exp(-Td*s)                    
           (1+Tp1*s)(1+Tp2*s)                                
                                                             
         Kp = 0.56884                                        
        Tp1 = 17989                                          
        Tp2 = 593.88                                         
         Td = 629.61                                         
                                                             
  From input "P_{pv}" to output "LOD":                       
                  Kp                                         
  G12(s) = ----------------- * exp(-Td*s)                    
           (1+Tp1*s)(1+Tp2*s)                                
                                                             
         Kp = 52.797                                         
        Tp1 = 37.875                                         
        Tp2 = 43828                                          
         Td = 8.79                                           
                                                             
  An additive ARMA disturbance model exists for output "LOD":
      y = G u + (C/D)e                                       
                                                             
      C(s) = s + 0.03899                                     
      D(s) = s + 0.002102                                    
                                                             
Name: T3s_2d
Parameterization:
    'P2D'    'P2D'
   Number of free coefficients: 10
   Use "getpvec", "getcov" for parameters and their uncertainties.

Status:                                                        
Estimated using PROCEST on time domain data "T3s_edit".        
Fit to estimation data: [91.07 94.18 92.84]% (prediction focus)
FPE: 0.1222, MSE: [0.1101 0.04204 0.07042]  

Is there a way to convert idproc to idss in canonical observable form?

[Rajiv Singh]

OK, thanks. Call FINDSTATES with Inf horizon, as in

x0 = findstates(sys, data, Inf)

since you are after simulation worthiness.

[julian navas]

That fixes it. Thank you very much again for your valuable help.