Error in simscape when computing derivative

I have a Simscape block that takes linPosition (physical signal) as an input and computes rotPosition (another physical signal). If I put it in my model, where linPosition is 0 for 1ms and then starts increasing smoothly, it works (rotPosition and theta are constant for 1ms and then follow the growth of linPosition). The problem arises when I try to add an additional computation: omega == der(theta) (I need it to make more calculation and provide another output). With that additional equation, I get this warning:

Equations (including nonlinear equations) of one or more components may be dependent or inconsistent. This can cause problems in transient initialization.

referred to the first and second equations. Then the simulation stops at 1ms:

An error occurred during simulation and the simulation was stopped

Caused by:

['testCircuit/Solver Configuration']: Transient initialization at time 0.001, solving for consistent states and modes, failed to converge.

Nonlinear solver: Linear Algebra error. Failed to solve using iteration matrix.

all components and nodal across variables involved

Cannot solve for one or more variables, including dynamic variable derivatives:

Time derivative of 'CRIB.CRIB_velocityBased.breaking_part.offsetSliderCrank.theta' (rot angle with respect to horizontal)

'CRIB.CRIB_velocityBased.breaking_part.offsetSliderCrank.omega' (omega)

In my mind, theta and gamma in the equations are computed and then omega follows the value of theta, but I guess that this is not how it works.

component offsetSliderCrank

inputs
    linPosition = {0, 'mm'};
end

outputs
    rotPosition = {0, 'deg'};
end

parameters
    lRod = { 45, 'mm' }; %
    rJunction =  { 10, 'mm' }; %
    xC = { 2, 'mm' }; %
    yC = { 40, 'mm' }; %

end


variables(Access=Protected)
    theta = {value={0,'deg'},imin={-90,'deg'},imax={89,'deg'}}; % rot angle with respect to horizontal
    gamma = {value={0,'deg'},imin={0,'deg'},imax={180,'deg'}}; % angle between lRod and horizontal
    omega = {0,'deg/s'};

end

equations
    lRod*cos(gamma)+xC == rJunction*cos(theta);
    yC-linPosition - rJunction*sin(theta) == lRod*sin(gamma);
    let
        theta0 = asin((yC^2+rJunction^2-lRod^2+xC^2)/(2*rJunction*sqrt(xC^2+yC^2))) - atan(xC/yC) ;
    in
        rotPosition == - (theta- theta0);
    end
    omega == der(theta);
end

end

Answer 1

[Torsten]

Can you compute der(linPosition) ?

If this is the case, you get an algebraic equation for omega = der(theta):

Write the first two equations as

lRod*cos(gamma) = rJunction*cos(theta) - xC
lRod*sin(gamma) = yC - linPosition(t) - rJunction*sin(theta)

Now square both equations and add them:

lRod^2 = (rJunction*cos(theta) - xC)^2 + (yC - linPosition(t) - rJunction*sin(theta))^2
lRod^2 =  -2*xC*rJunction*cos(theta) + xC^2 + (yC - linPosition(t))^2 - 2*(yC - linPosition(t))*rJunction*sin(theta) + rJunction^2
lRod^2 - xC^2 - yC^2 - rJunction^2 = -2*xC*rJunction*cos(theta) - 2*yC*linPosition(t) + linPosition(t)^2 - 2*(yC - linPosition(t))*rJunction*sin(theta)

Now differentiate both sides with respect to t:

0 = 2*xC*rJunction*sin(theta)*der(theta) - 2*yC*der(linPosition) + 2*linPosition*der(linPosition) - 2*yC*rJunction*cos(theta)*der(theta) + 2*linPosition*rJunction*cos(theta)*der(theta) + 2*der(linPosition)*rJunction*sin(theta)

and solve for der(theta):

der(theta) = (2*yC*der(linPosition) - 2*linPosition*der(linPosition) - 2*der(linPosition)*rJunction*sin(theta))/(2*xC*rJunction*sin(theta)+2*yC*rJunction*cos(theta)+2*linPosition*rJunction*cos(theta))

By the way: if gamma and theta are in degrees, is it possible in Simscape to use sin and cos ? In MATLAB, you had to use sind and cosd.

[Ilario Triscari]

Thank you, the workaround works well, but still I don't understand why omega == der(theta) doesn't

[Torsten]

Thank you, the workaround works well, but still I don't understand why omega == der(theta) doesn't

Most probably because "theta" is not given explicitly, but implicitly by a system of algebraic equations.

To test this idea, try whether it works if you solve

lRod^2 - xC^2 - yC^2 - rJunction^2 = -2*xC*rJunction*cos(theta) - 2*yC*linPosition(t) + linPosition(t)^2 - 2*(yC - linPosition(t))*rJunction*sin(theta)

explicitly for theta (by replacing sin(theta) by +/- sqrt(1-cos^2(theta)), isolating +/- sqrt(1-cos^2(theta)) on one side of the equation, squaring both sides and solving the resulting quadratic equation in cos(theta) for cos(theta)) and then apply omega = der(theta).

In case you define "theta" explicitly, you will have to skip one of the two equations defining "theta" and "gamma":

lRod*cos(gamma)+xC == rJunction*cos(theta);
yC-linPosition - rJunction*sin(theta) == lRod*sin(gamma);

[Torsten]

The following simple example shows that your system of differential-algebraic equation will be hard to solve in its original form.

The problem is that

der(theta) = omega

is interpreted as a differential equation for "theta", not as an evaluation of the time derivative of "theta" that can explicitly be computed from the former two algebraic equations for "theta" and "gamma". Thus the solver assumes two equations for "theta" and no equation for omega: one algebraic and one differential equation. This results in an index problem as error message (interprete y(1) as "theta", y(2) as "gamma" and y(3) as "omega"):

fun = @(t,y)[y(1)-y(2)+1-sin(t);3*y(1)+4*y(2)-7;y(3)];
mass = [0,0,0;0,0,0;1,0,0];
tspan = [0 1];
y01 = [1 -1;3 4]\[-1;7];
y02 = 4/7;
y0 = [y01;y02];
options = odeset('Mass',mass);
[T,Y] = ode15s(fun,tspan,y0,options)

Error using daeic12 (line 72)
This DAE appears to be of index greater than 1.

Error in ode15s (line 204)
[y,yp,f0,dfdy,nFE,nPD,Jfac] = daeic12(ode,odeArgs,t,ICtype,Mt,y,yp0,f0,...

plot(T,Y)

So even if you explicitly solve for "theta" as suggested in my above response, I doubt that you will overcome the problem because the structure of your system of differential-algebraic equations remains the same:

theta = f1(gamma)
f2(theta,gamma) = 0
dtheta/dt = omega