My outputs are unrealistically high and I have no idea why, d,x,xdot,xdotdot at in magnitudes beyond what it needs to be

Delta_t=1;
theta=1.4;
a_0=6/(theta*Delta_t)^2;
a_1=3/(theta*Delta_t);
a_2=a_1*2;
a_3=(theta*Delta_t)/2;
a_4=a_0/theta;
a_5=-a_2/theta;
a_6=1-3/theta;
a_7=Delta_t/2;
a_8=(Delta_t^2)/6;

%Stiffness Matrix

M=22; %lb/ft
K=0.8; %Stiffness Coeff.
C=0.25674; %Damping Coeff.

K_hat=a_0*M+a_1*C+K

K_hat = 68.6971

%Payload

L=100;
t=1:300;
Mw = zeros(numel(t),1);
Mw(t<80) = 20;
Mw(t>=80 & t<=90) = 35;
Mw(t>90 & t<=100) = 60;
Mw(t>100 & t<=200)= 100;
Mw(t>200) = 20;
Fv=0.04*L*Mw; %Variable of the Load
F=Fv.';

x=3; %ft
x_dot=3; %ft/s
x_dotdot=3; %ft/s^2
for ii=1:220
    F_hat(ii+1)=F(ii)+theta*(F(ii+1)-F(ii))+M*(a_0*x(ii)+a_2*x_dot(ii)+2*x_dotdot(ii))+C*(a_1*x(ii)+2*x_dot(ii)+a_3*x_dotdot(ii));
    d(:,ii+1)=F_hat(:,ii)/K_hat; % x=F/K
    x_dotdot(:,ii+1)=a_4*(d(:,ii+1)-x(:,ii))+a_5*x_dot(:,ii)+a_6*x_dotdot(:,ii);
    x_dot(:,ii+1)=x_dot(:,ii)+a_7*(x_dotdot(:,ii+1)+x_dotdot(:,ii));
    x(:,ii+1)=x(:,ii)+Delta_t*x_dot(:,ii)+a_8*(x_dotdot(:,ii+1)+2*x_dotdot(:,ii));
end
plot(x_dot)

[Walter Roberson]

How does this differ from https://www.mathworks.com/matlabcentral/answers/1670259-i-am-unable-to-get-an-array-of-x-x_dot-and-x_dotdot-for-a-series-of-different-f-inputs-down-at-the?s_tid=srchtitle ?

[Mark Johnson]

It doesnt except I asked a different question. One being I couldnt get it to work and the second being its working but now with unrealistic outputs.

[Mark Johnson]

Yeah the outputs become insanely high and I want to be able to keep them stable

[Walter Roberson]

Which of the two questions do you want answered? You should close the other one.

[Mark Johnson]

I would like to figure out why the outputs are so unrealistically huge. I will close out the other question ticket.

[the cyclist]

(Edited question so that the 220 iterations worth of output are not displayed.)

Answer 1

[Walter Roberson]

Your outputs wobble a lot; they just look straight because the wobbles increase enough as you go that the previous part is straight by comparison .

Switching to symbolic shows that the problem is not with round-off error in the calculations: there is something in the formulas that leads to this state.

Q = @(v) sym(v);
Delta_t = Q(1);
theta = Q(1.4);
a_0=6/(theta*Delta_t)^2;
a_1=3/(theta*Delta_t);
a_2=a_1*2;
a_3=(theta*Delta_t)/2;
a_4=a_0/theta;
a_5=-a_2/theta;
a_6=1-3/theta;
a_7=Delta_t/2;
a_8=(Delta_t^2)/6;

%Stiffness Matrix

M = Q(22); %lb/ft
K = Q(0.8); %Stiffness Coeff.
C = Q(0.25674); %Damping Coeff.

K_hat=a_0*M+a_1*C+K

K_hat = 

%Payload

L = Q(100);
t = Q(1:300);
Mw = zeros(numel(t),1,'sym');
Mw(t<80) = 20;
Mw(t>=80 & t<=90) = 35;
Mw(t>90 & t<=100) = 60;
Mw(t>100 & t<=200)= 100;
Mw(t>200) = 20;
Fv = Q(0.04)*L*Mw; %Variable of the Load
F=Fv.';

x = Q(3); %ft
x_dot = Q(3); %ft/s
x_dotdot = Q(3); %ft/s^2
for ii=1:220
    F_hat(ii+1) = F(ii)+theta*(F(ii+1)-F(ii))+M*(a_0*x(ii)+a_2*x_dot(ii)+2*x_dotdot(ii))+C*(a_1*x(ii)+2*x_dot(ii)+a_3*x_dotdot(ii));
    d(ii+1) = F_hat(ii)/K_hat; % x=F/K
    x_dotdot(ii+1) = a_4*(d(ii+1)-x(ii))+a_5*x_dot(ii)+a_6*x_dotdot(ii);
    x_dot(ii+1) = x_dot(ii)+a_7*(x_dotdot(ii+1)+x_dotdot(ii));
    x(ii+1) = x(ii)+Delta_t*x_dot(ii)+a_8*(x_dotdot(ii+1)+2*x_dotdot(ii));
end
plot(double(x_dot))

x_dot(end)

ans = 

vpa(ans)

ans = 

1.0991241567310856704335408598327e+88

plot(double(x_dot(1:150)))

plot(double(x_dot(1:50)))

[Mark Johnson]

So I will need to add a stabilization matrix of some form?

[Walter Roberson]

Sorry, I do not know anything about stabilization matrices.

[Torsten]

What integration method do you use for your equations ? I can't figure it out from your updates in the loop.

[Mark Johnson]

These are the equations I am trying to make code for and run show that it shows the dynamic plots

[Torsten]

I don't know if you can deduce the updating equations from above.

If yes, you should do that because the risk of typos is always quite high.

Assuming the equations are correct, the three equations under (2.3) are implicit in xdotdot(t+theta*dt), xdot(t+theta*dt) and x(t+theta*dt). So at least you will have to solve a linear system to solve for them (which I cannot see in your code).

[mark johnson]

This is the linear Equations I was able to get but I don’t know how to implement into matlab

[Mark Johnson]

[Torsten]

And because you don't know how to implement something you just program something else ?

At least you could have told us already in your question.

x(1) = 3;
xdot(1) = 3;
xdotdot(1) = 3;
for i=1:220
    A = [1 0 -a_4; -a_7 1 0; -a_8 0 1];
    b = [-a_4*x(i)+a_5*xdot(i)+a_6*xdotdot(i);xdot(i)+a_7*xdotdot(i);x(i)+Delta_t*xdot(i)+2*a_8*xdotdot(i)];
    sol = A\b;
    xdotdot(i+1) = sol(1);
    xdot(i+1) = sol(2);
    x(i+1) = sol(3);
end

Whether this is correct - I don't know. I just took the expressions from your sheet.

[Mark Johnson]

I have tried numerous tactics to make the research documents equations work with a foor loop. Linearization is just one thing I considered before going with the set up shown before. I understand matrix b but what is matrix A? And can this be implemented with the changing F_hat?

[Torsten]

I wrote the equations shown as

A*[xdotdot(i+1);xdot(i+1);x(i+1)] = b

where A is a 3x3 matrix and b is a 3x1 vector and solved for [xdotdot(i+1);xdot(i+1);x(i+1)].

From your sheet, I could not deduce how F_hat comes into play here. Equations (2.3) are not dependent on F_hat, as far as I can see.

But you are the person who has the background to change A and b adequately.

[Mark Johnson]

The relation to F_hat comes in the equation x(ii+1)=F_hat/K_hat where K_hat is just a value. From the image I sent u(k) of the state space equation is [F_hat/K_hat]. This means that the three equations are dependent on F_hat.

[Sam Chak]

Hi @Mark Johnson

When things like this (instability or singularity) happen, most of the time it is because

1. the system is naturally unstable;
2. the system is modeled incorrectly, either failing to capture some critical dynamics, or using the equation that was errorly printed on the reference material (usually only experts can spot this);
3. the structure of the system is basically correct, but due to oversight, some ± signs in the terms are switched unintentionally (i.e., the dotdot term stays on the left, and the rest goes to the right-hand side).

For a typical rigid-flexible coupling system like the , most likely it is Case #3. Here are two examples:

Stable:

Unstable:

and it takes only one sign to destabilize the entire system.

What we can help:

You show us the model or equations from some reputable references, and we can check.

What we cannot help (with tips):

If you derive the model from a set of fundamental laws, then you need to check it yourself, unless there is an expert who can check whether your derivations steps are correct or not. Else, try using a different modeling approach (say, formulate the model using the Lagrangian mechanics), and see if you arrive at the same model. If it is the same, then you can probably argue that this is Case #1 (naturally unstable).

[Mark Johnson]

P_2492.pdf (sensorsportal.com) Here is the link of the research study done. This is what I am basing this code on. I have made progress. Its now stable and corrects itself.

[Sam Chak]

Hi @Mark Johnson, it's to good to hear that you've made a good effort to check and correct it. 👍