Solve closed loop system equation with disturbance input via function handle

Hi all,

i have the following issue: I want to simulate a closed loop system behaviour and therefore I want to solve the ode with a function handle. The closed loop experiences a disturbance input which is a matrix where rows corresoponds to disturbance inputs, and columns to the time steps of the simulation. This means the columns of the disturbance matrix are equal to the time steps of the simulation. How can I simulate this, especially without a for loop where I dont have to specify for each time step the correct disturbance input column?

Here is some code:

[~,xCL] = ode15s(@(t,x,z)closedLoop(t,x,zFull,sysSSFULL),t_span,zeros(592,1));
%zFull is full disturbance matrix with disturbance vector for all timesteps, zeros(592,1) is just initial condition

corresponding function is

function dx = closedLoop(t,x,z,sys)

dx = (sys.A - sys.B_u*sys.Kopt)*x+z ; ... %z is disturbance input vector at corresponding time
    
end

I always get either the error that matrices are not consistent, which would mean matlab takes the whole matrix instead of the correct vector for each time step, or not enough input arguments if I change the above declaration around when calling the anonymous function. What am I making wrong?

Thanks for helping!

[Sam Chak]

Hi @Fabian Friz,

The initial condition syntax 'zeros(592, 1)' indicates that there are 592 state variables. This seems like a very large system. Also, the example in my answer assumes matched disturbances, where the external disturbance terms appear at the same level of differentiation as the system input.

Answer 1

[Walter Roberson]

In the case where the disturbances are considered instantaneous in duration, or instantaneous in starting up, then the answer is that you do not do that with any of the ode*() functions -- not unless you terminate the ode*() call right when the disturbance would start and then call the ode*() function again (with possibly adjusted boundary conditions) . The mathematics of Runge-Kutta methods is only valid when the function you provide has continuous second derivatives.

Using interp1() with the default 'linear' option does not satisfy the mathematics. Using interp1() with a cubic spline does satisfy the boundary conditions.

However, if you are doing interpolation inside the function then it logically must be a situation where you would be fine with the code "looking forward in time" to see what the disturbance is going to be, and already being in motion towards that value in the steps near the disturbance. If the situation at hand is, for example, irregular rocks in the road, then the car springs cannot realistically be moving to react to a rock that has not been reached yet. However, if the situation at hand were a series of speed-bumps of different height but consistent curvature, then by knowing the location and the curvature you could project back to the point at which the tire will first encounter the smooth edge of the bump, and from that location to the center of the bump you could climb smoothly, moving smoothly to reach the peak -- and in such a case a good interpolation function might logically work.

Answer 2

[Torsten]

By using interp1 to interpolate the column value to the time instant of the ODE integrator.

Or use your own integrator with a fixed time step right from the beginning.

Answer 3

[Sam Chak]

Hi @Fabian Friz

In addition to using the 1-D data interpolation method interp1(), as suggested by @Torsten, it is also possible to employ a curve-fitted model if the goodness-of-fit of the disturbance model is satisfactory within the simulation time interval. Here is an example:

% time-dependent disturbance, d(t)
t      = 1/3*[0:5:90]';
d      = [0 0.1 sqrt(3)/10 0.2 sqrt(3)/10 0.1 0 -0.1 -sqrt(3)/10 -0.2 -sqrt(3)/10 -0.1 0 0.1 sqrt(3)/10 0.2 sqrt(3)/10 0.1 0]';
model  = fittype('a*sin(pi/b*t)', 'dependent', {'d'}, 'independent', {'t'}, 'coefficients', {'a', 'b'});
[myfit, gof] = fit(t, d, model, 'start', [0.25, 12.5])

myfit = 
     General model:
     myfit(t) = a*sin(pi/b*t)
     Coefficients (with 95% confidence bounds):
       a =         0.2  (0.2, 0.2)
       b =          10  (10, 10)

gof = struct with fields:
           sse: 1.1023e-12
       rsquare: 1.0000
           dfe: 17
    adjrsquare: 1.0000
          rmse: 2.5464e-07

plot(myfit, t, d), grid on

% ODE solver
tspan  = [0 30];
x0     = [1  0];
[t, x] = ode15s(@(t,x,z) closedLoop(t, x, myfit), tspan, x0);
plot(t, x), grid on, xlabel('Time'), legend('x_{1}', 'x_{2}')

% dynamic system
function dx = closedLoop(t, x, myfit)
    A    = [0  1; ...   % state matrix
            0 -1];
    B    = [0; ...      % input matrix
            1];
    Kopt = [1  1];      % optimal feedback gain matrix
    d    = myfit.a*sin(pi/myfit.b*t);   % fitted disturbance model
    z    = [0; ...      % disturbance vector
            d];
    dx   = (A - B*Kopt)*x + z;
end