How can I integrate state-space equations in 6DOF with ode45?

I have a state-space representation of a system in six degrees of freedom, and have used the forward Euler method to update the state, as seen in the attached script. However, this gave an inaccurate result, so I want to change to the Runge-Kutta method. I have tried to use the ode45 function, but am unsure about the syntax.

Here is my attempt (with the rest of the code unchanged):

    [t,x]=ode45(@(t,x) f(t, x, N, u, A, B),t,x0);

    for k=1:N-1 % need -1 here so that x(n+1) does not lead to an index>N
        for i = 1:1:12
            y(i,k+1)=C(i,:)*x(:,k+1)+D(i,:)*u(:,k+1);       % y is an output based on x and u, so is not integrated
        end
    % Calculate new version of A
    A=calcA(x(:,k+1));
    end

    function [t, x]=f(t, x, N, u, A, B)
    for k=1:N-1 
        for i = 1:1:12
            x(i,k+1)=A(i,:)*x(:,k)+B(i,:)*u(:,k);
        end
    end
  end

I get the following errors:

Error using odearguments (line 95) @(T,X)F(T,N,X,U,A,B) returns a vector of length 1, but the length of initial conditions vector is 12. The vector returned by @(T,X)F(T,N,X,U,A,B) and the initial conditions vector must have the same number of elements.

Error in ode45 (line 115) odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

Error in newsys4 (line 63) [t,x]=ode45(@(t,x) f(t, x, N, u, A, B),t,x0);

I think the correct length is 12, since there are 12 states. How can I make @(T,X)F(T,N,X,U,A,B) return a vector of length 12?

[Jan]

The function, which is integrated by ODE45 replies one output only: the derivative as vector. dx = f(t, x, ...) . It is not clear to me, which operations you apply and why a matrix is returned.

[EllaM]

The ode45 function should return a vector with 12 elements (since there are 12 states), but apparently it returns a vector with only 1 element, and I'm not sure why.

Answer 1

[James Tursa]

Your derivative function f should only be returning the derivative of x, not t also. I.e., the signature of f should be something like this:

    function dx = f(t, x, N, u, A, B)
    dx = _______; % Add code here to calculate the derivative of x at time t
    end

The way you have it now, ode45 is simply getting as a result of the f call the first returned argument, t, which is a scalar (1 element). Hence the error.

[EllaM]

Sorry, but I still can't make it work. I have changed the function and function call like this:

    [x]=ode45(@(x,u) f(t, x, N, u, A, B),t,x0);

    function [x]=f(t, x, N, u, A, B)
        for k=1:N-1 
            for i = 1:1:12
                x(i,k+1)=A(i,:)*x(:,k)+B(i,:)*u(:,k);
            end
        end
    end

This gives me the following errors:

Error using * Inner matrix dimensions must agree.

Error in f (line 5) x(i,k+1)=A(i,:)*x(:,k)+B(i,:)*u(:,k);

Error in newsys4>@(x,u)f(t,x,N,u,A,B)

Error in odearguments (line 90) f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.

Error in ode45 (line 115) odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

Error in newsys4 (line 63) [x]=ode45(@(x,u) f(t, x, N, u, A, B),t,x0);

I am not sure why the matrix dimensions don't agree, since the dimensions are as follows:

A:12x12, x:12xN, B:12x6, u:6xN

Thanks!

[Jan]

@EllaM: Do you see, that James' suggestion replies dx and not the input x?

Use the debugger to find the dimensions:

dbstop if error

Now run the code until it stops at the error. Now check what the dimensions are, not what you expect them to be.

[EllaM]

I see, so the problem is that x and u are not passed correctly to the function f. Is there something wrong in the function call?

    [xdot]=ode45(@(x,u) f(t, x, N, u, A, B),t,x0);

    for k=1:N-1 
        for i = 1:1:12
            x(i,k+1)=xdot(i,k);
            y(i,k+1)=C(i,:)*x(:,k+1)+D(i,:)*u(:,k+1);
        end

        A=calcA(x(:,k+1));
    end

    function [xdot]=f(t, x, N, u, A, B)
        xdot = zeros(12,N);
        for k=1:N 
            for i = 1:1:12
                xdot(i,k)=A(i,:)*x(:,k)+B(i,:)*u(:,k);
            end
        end
    end

[Jan]

xdot = ode45(@(t,x) f(t, x, N, u, A, B), t, x0)

not @(x,u). And xdot must be a vector, if you reply it to ode45. It must have the same dimensions as the inital value x0, because it is the derivative ate the time t and position x.

[EllaM]

Ah, thank you, that makes sense!