Best way to solve differential algebraic equations?

Hello, I went through the pendulum example about solving differential algebraic equations. The system in the example is an autonomous system. For my application I have to solve an actuated system, so for the purpose of learning, I tried to modify the pendulum system so that it is actuated. To make it simple, I just made the pole radius r time variable (r -> r(t) = 1/10*(cos(t)+1) + 3). This is in line 44 in my code. The problem is, that decic gives me an error "Index exceeds matrix dimensions."

Can someone tell me how I can include a time dependent actuation function in a way, so I can change the function easily? Preferrably, I only want to have to change r(t) = 1/10*(cos(t)+1) + 3 to something else without having to change anything else in the code.

Thank you and regards.

[Torsten]

Could you include the code you are using as an ASCII file ?

[Matthias Scheuerer]

Hi Torsten, I attachted the code in my answer below.

[Torsten]

clear all; close all;

syms x(t) y(t) T(t) m r(t) g;

% system of equations
eqn1 = m*diff(x,t,2) == T*x/r;
eqn2 = m*diff(y,t,2) == T*y/r - m*g;
eqn3 = x^2+y^2 == r^2;
eqn4 = r == 1/10*(cos(t)+1) + 3;
eqns = [eqn1 eqn2 eqn3 eqn4]

% vector of state variables
vars = [x(t); y(t); T(t); r(t)]

% dimension of unmanipulated state vector
origVars = length(vars)

M = incidenceMatrix(eqns, vars)

% reduce to first order system
[eqns, vars, R] = reduceDifferentialOrder(eqns, vars)

incidenceMatrix(eqns, vars)

% check if DAE is of order 0 or 1
isLowIndexDAE(eqns, vars)

% reduce index to 0 or 1
[DAEs, DAEvars, DAER] = reduceDAEIndex(eqns, vars)
incidenceMatrix(DAEs, DAEvars)

% eliminate redundant state variables
[DAEs, DAEvars, DAER] = reduceRedundancies(DAEs, DAEvars)
incidenceMatrix(DAEs, DAEvars)

% now system has index 0 or 1
isLowIndexDAE(DAEs,DAEvars)

% determine parameter vector
pDAEs = symvar(DAEs)
pDAEvars = symvar(DAEvars)
extraParams = setdiff(pDAEs, pDAEvars)
DAEs = subs(DAEs)
% create function handle for solver
f = daeFunction(DAEs, DAEvars, g, m)

% set numerical values for parameters
g = 9.81;
m = 1;
%r = 1;

% create function handle only containing numerical parameters for solver
F = @(t, Y, YP) f(t, Y, YP, g, m)

DAEvars

My guess is that in the sequel, you will have to supply 8 instead of 7 initial conditions.

But just test it.

Best wishes

Torsten.

[Matthias Scheuerer]

Hi Torsten, I tried your suggestion, but the error is still the same. Actually the state vector (after all the previous modifications) is

DAEvars = [x(t); y(t); T(t); Dxt(t); Dyt(t); Dytt(t); Dxtt(t)],

so it contains 7 elements, hence 7 intial values should be enough, right? The system of equations also contains 7 equations.

[Torsten]

The "dirty" method is to use the old code, remove r from the syms variables and replace "r" in equations eqn1 - eqn3 by "1/10*(cos(t)+1) + 3".

Answer 1

[Matthias Scheuerer]

Matlab Support told me the solution: The problem is, that r(t) is a symbolic function, so r(0) returns a symbolic number. I used r(0) in the expression for a guess for inital values. This causes decic to throw an error, because it can only handle numerical inputs. So I had to cast y0est into double (double(y0est)) and it works. However, "Index exceeds matrix dimensions." is not a very suitable error message in this case.

Answer 2

[Matthias Scheuerer]

Hi Birdman, thank you for your suggestion. r(t) is already a symbolic expression in my script. I subsitute it with the expression above, so there is t as a symbolic variable left. But it still produces the error I described above.

I attached the code as ASCII file.