1

For training purpose, I try to modeling and analyze a dynamic system using python. I follow the example given on this site.

In this example, they are using ODEint to resolve the differential equation and simulate a step response. This method gives the following result: enter image description here

Then, I'm using the state space representation to simulate the same system. This gives the following code:


import scipy.signal as signal
import numpy as np
import matplotlib.pyplot as plt

m1 = 1.2
m2 = 1.2
k1 = 0.9
k2 = 0.4
b1 = 0.8
b2 = 0.4

A = [[0, 0, 1, 0],[0, 0, 0, 1],[(-k1-k2)/m1, k2/m1, -b1/m1, 0],[k2/m2, -k2/m2, 0, -b2/m2]]
B = [[0],[0],[1],[0]]
C = [[1, 0, 0, 0]]
D = [0]

sys = signal.StateSpace(A, B, C, D)
x0 = [0,0,0,0]
start = 0
stop = 40
step = 0.01
t = np.arange(start, stop, step)
u = np.ones(len(t))
t, y, x = signal.lsim(sys, u, t, x0)

plt.plot(t, x)
plt.legend(['m1_position', 'm2_position', 'm1 speed', 'm2_speed'])
plt.grid(True)

And the result of the state space simulation: enter image description here

We see that the results are different between the two methods. The shapes are identical, but the value are different. I expect that the both methods give the same values. I have tried to understand if the both methods use different resolution algorithm, but lsim documentation do not give any details for that.

I would like to understand the reason for this difference? Is one method more accurate than the other? In this particular case or in general? How can we improve the accuracy?

Max
  • 11
  • 3

1 Answers1

1

The equilibrium position in your coded system is at x2 = x1 = m1/k1 = 1.2/0.9=4/3=1.3333, which is where the solution goes to.

To get the same behavior as the explicitly coded ODE function, you need to set B = [[0],[0],[1/m1],[0]] to correctly translate force into acceleration, to then get the equilibrium at x2 = x1 = 1/k1 = 1/0.9=1.111.

As the force input is the only non-homogeneity, the solutions differ, as observed, only by a scale factor.

Lutz Lehmann
  • 25,219
  • 2
  • 22
  • 51
  • Oh yes ! So the point is that I made a mistake writing the B matrix. B matrix is not `B = [[0],[0],[1],[0]]` but `B = [[0],[0],[1/m1],[0]]`. Thank you ! – Max Jan 23 '23 at 23:03