# How do I stabilize my oscillating system?

I have two scalars $x$ and $y$ that vary with time $t$ such that

\begin{align} \frac{dx}{dt} &= ay + b + dx\\ \frac{dy}{dt} &= cx \end{align} but $a, b, c$ are unknown.

If $d$ is too small, the system will be underdamped (and oscillate badly), but if $d$ is too large, the system is overdamped, and doesn't converge fast enough.

How do I set $d$ based on statistics taken about historic $x, y$ so that the system is critically damped and moves as quickly as possible to the point where $x=0, y=\bar y = -\frac{b}{a}$?

• Note that your system cannot be critically damped because you have no damping at all. Your differential equation is $y''-acy=bc$. For damping you would need a term $y'$. What you can do is play around with the initial conditions to get rid of the unstable pole (by pole-zero cancellation) which results in an exponential decay with time constant $1/\sqrt{ac}$ (assuming $ac>0$). That's what is explained in Ced's answer. – Matt L. Jan 31 '18 at 15:35
• @MattL. My question is actually: how do I solve for that $y'$ term that critically damps my system? I'll comment on the answer too. – Neil G Jan 31 '18 at 16:48
• Your differential equation would need to be of the form $y''+ay'+b=c$ (all constants unequal zero) in order to make critical damping possible, but it isn't, that was my point. – Matt L. Jan 31 '18 at 17:13
• @Niel G, including a $d \cdot y$ term in your second equation would do the trick. The coefficient $d$ would have to have a negative value or your system would grow instead of decay. (Negative feedback loop). Since you found oscillation in your gradient solution, I suspect your $ac < 0$. The nature of your solution has to do with the roots of $r^2 + Ar + B = 0$ from Matt L's differential equation. I capitalized the coefficients to distinguish them from yours. – Cedron Dawg Jan 31 '18 at 18:30
• @Niel G. (continued) If $ac > 0$, then the solution I gave in my answer is unstable. The initial conditions have to be met exactly or C1 will have a non-zero value and your system will go out of control. – Cedron Dawg Jan 31 '18 at 18:30

First, combine your two variable set of first order differential equations into a single variable second order one.

$$\frac{d^2y}{dt^2} = c \frac{dx}{dt} = acy + bc$$

$$\frac{d^2y}{dt^2} - acy = bc$$

Solve:

$$y(t) = C_1 e^{\sqrt{ac} \cdot t} + C_2 e^{-\sqrt{ac} \cdot t} - \frac{b}{a}$$

$$x(t) = \int{ ( ay(t) + b ) dt }$$

$$x(t) = C_1 \frac{a}{\sqrt{ac}} e^{\sqrt{ac} \cdot t} - C_2 \frac{a}{\sqrt{ac}} e^{-\sqrt{ac} \cdot t} + C_3$$

Now analyze:

Initial conditions have to be set so that $C_1 = 0$ or your system will grow exponentially (assuming $ac > 0$). Also, $C_3$ will need to be zero. If $ac < 0$ then the system will oscillate forever since there is no damping term in the second order equation.

$$\frac{dy}{dt} = \sqrt{ac} \left( C_1 e^{\sqrt{ac} \cdot t} - C_2 e^{-\sqrt{ac} \cdot t} \right)$$

$$\frac{dx}{dt} = a \left( C_1 e^{\sqrt{ac} \cdot t} + C_2 e^{-\sqrt{ac} \cdot t} \right)$$

Initial Conditions:

$$y(0) = C_1 + C_2 - \frac{b}{a}$$

$$\frac{dy}{dt}(0) = \sqrt{ac} \left( C_1 - C_2 \right)$$

$$x(0) = \frac{a}{\sqrt{ac}} ( C_1 - C_2 ) + C_3$$

$$\frac{dx}{dt}(0) = a \left( C_1 + C_2 \right)$$

With $C_1, C_3 = 0$:

$$y(0) = C_2 - \frac{b}{a}$$

$$\frac{dy}{dt}(0) = - \sqrt{ac} \cdot C_2$$

$$x(0) = -\frac{a}{\sqrt{ac}} \cdot C_2$$

$$\frac{dx}{dt}(0) = a \cdot C_2$$

It's been a while, I hope I did it right.

Ced

==================================

Followup:

Yes, you can include an $x$ term in the first equation, or a $y$ term in the second equation, or both. Then repeat the process I went through above to put your system in this form:

$$\ddot y + A \dot y + B y = 0$$

Whatever constant you get on the right side will be part of the particular solution. Add it in later as I did above.

This is the basics of solving a homogenous second degree linear differential equation:

Assume the solution is of the form $C e^{rt}$

$$r^2 C e^{rt} + A r C e^{rt} + B C e^{rt} = 0$$

$$\left( r^2 + A r + B \right) \left( C e^{rt} \right) = 0$$

One, or both, of the factors has to be zero. $C e^{rt}$ won't be, so:

$$r^2 + A r + B = 0$$

This is a standard quadratic equation. With real valued coefficients you will either get two real roots, or a conjugate pair of complex roots.

$$r = - \frac{ A }{2} \pm \sqrt{ \left(\frac{ A }{2} \right)^2 - B }$$

This will give you $r_1$ and $r_2$. Since it is a linear differential equation, any linear combination of the two solutions will also be a solution.

$$y(t) = C_1 e^{r_1 t} + C_2 e^{r_2 t}$$

The real part of the root will determine the growth behavior. If it is negative the system will converge to an equilibrium point, if it is positive, the system will grow exponentially. The imaginary part, if there, will determine the frequency of the oscillations. Therefore, you can calculate the parameters from observations by measuring the rate of decay to get the real part(s), and measuring the frequency to get the imaginary part. Once you have those, you should be able to solve for your $a,b,c$ and $d,e$ if you include them.

Of course, you could leave your equation set in matrix form and solve it with eigenvalues and eigenvectors, but that is a whole 'nuther discussion.

Ced

• Makes sense, thanks. This looks like the "damped oscillator" I've been recently reading about. I've added a comment to explain what I'm really after. – Neil G Jan 31 '18 at 23:52
• I understand everything you added. I wish I could +1 again. I've tried to clarify exactly what I'm after. Sorry for my roundabout way of asking, but I am still learning this. – Neil G Feb 1 '18 at 0:52

An empirical approach:

I'm adding a new answer as it is completely different in nature.

It seems you have a system you are studying where you can vary one parameter (or more?) to control the dampening rate, and you are looking for the best setting, i.e. finding the critical dampening. This is where the roots of the second order differential are the same, the boundary between oscillating solutions and decay solutions.

Simply measure the frequency of your system under several settings of your dampening factor, building a table like this:

$$\begin{array}{cc} Parameter & Frequency \\ p_1 & f_1 \\ p_2 & f_2 \\ p_3 & f_3 \\ . & . \\ . & . \\ . & . \\ p_N & f_N \\ \end{array}$$

Then plot the values on a graph. Extrapolate the curve (manually or algorithmically) to where it intersects the $p$ axis. That is the parameter value you want.

Ced

P.S. One wrinkle that may or not may have an impact, when the solution of the characteristic equation is a double root, the general solution of the equation changes to:

$$y(t) = C_1 e^{rt} + C_2 t e^{rt}$$

========================================

If you carry your new formulation into a second order differential equation of $y$ you get:

$$\frac{d^2y}{dt^2} - d \frac{dy}{dt} - acy = bc$$

The homogenous portion is:

$$\ddot{y} + A \dot{y} + B y = 0$$

So:

$$A = -d$$

$$B = -ac$$

The frequency of your oscillations (in cycles per unit time) will be:

$$f = \frac{1}{2\pi} \sqrt{ B - \left(\frac{ A }{2} \right)^2 }$$

As you vary your $d$ (and thus $A$) towards overdampening, the frequency will approach zero, reaching zero at the critical value. The graph is the square root of a convex downward parabola so there will be an intercept. [Edit: Doh, it will be the top half of an ellipse.] There isn't any asymptotic behavior as you ask in your first comment. It is the value of the function, not the slope that matters. Beyond that point, both roots are real and your solution is the sum of two decays, the frequency is imaginary.

Solving for $\frac{dy}{dA}$ can be done from the equations I gave in my previous answer, but I don't see how it will help find the critical point.

Ced

=======================================

Appendum:

You can use $x$ for the frequency measurement just as well as $y$, but you knew that.

The frequency won't change as the signal decays so you can measure the frequency as early as possible before the signal decays to where the noise is significant compared to the amplitude.

The good news is you really only have to calculate one calibration point, compared to the many I recommended above, since the expected equation is known.

$$4\pi^2 f^2 + \left(\frac{ A }{2} \right)^2 = B$$

Let $A_m$ be the parameter setting you are measuring the frequence for, and $f_m$ be the resulting frequency.

$$B = 4\pi^2 f_m^2 + \left(\frac{ A_m }{2} \right)^2$$

Setting $f=0$ to find your critical $A_c$:

$$\left(\frac{ A_c }{2} \right)^2 = 4\pi^2 f_m^2 + \left(\frac{ A_m }{2} \right)^2$$

Now solve for $A_c$:

$$A_c = 2 \cdot \sqrt{ 4\pi^2 f_m^2 + \left(\frac{ A_m }{2} \right)^2 }$$

You may want to take several measurements anyway to see if the results are consistent as expected.

• It makes sense that as $p \to 0$, the frequency goes to some number, say $f_0$, and as $p \to \infty$, the frequency goes to zero since an over-damped system does no oscillate. I don't see where there is any intercept? – Neil G Feb 1 '18 at 18:38
• I've also been thinking about this, and my feeling is that we should take the derivative of the velocity of $x$ (or $y$) wrt to the parameter $d$ (or $p$ as you've called it here). The velocity wrt to $d$ should be convex, so now we do have an intercept where the velocity does not change when $d$ changes, and that's where we want to stop. – Neil G Feb 1 '18 at 18:40
• Regarding your update, it's very hard afaik to measure the frequency of a signal. If the signal is noisy, then what will you say its frequency is? Isn't this doable with velocity instead? – Neil G Feb 1 '18 at 20:16
• @Neil G, In a situation like this you can smooth your signal and simply calculate the distance between zero crossings. More accuracy can be achieved by measuring across several oscillations. – Cedron Dawg Feb 1 '18 at 21:00
• I'm worried that if $x$ converges to zero, it's going to look like an infinite number of zero crossings, and then $d$ will go to $\infty$… – Neil G Feb 1 '18 at 21:04

Converting the model to a state space model

$$\begin{pmatrix}\dot x\\\dot y\end{pmatrix} = \begin{pmatrix}d &a\\ c &0\end{pmatrix}\begin{pmatrix}x\\y\end{pmatrix} + \begin{pmatrix}1\\ 0\end{pmatrix}b = A\hat x + Bb$$

Hence you would like to have the eigenvalues of the $A$ matrix to be on the open left half plane which is given already

$$\lambda^2 -d\lambda -ac = 0 \Longrightarrow \lambda_{1,2} = \frac{-d \pm \sqrt{d^2 + 4ac}}{2}$$ Obviously, these should first lie on the openleft half plane for stability. So that gives a first set of constraints on $d,a,c$.

For stable eigenvalues, if you further would like to match these to a typical second order system oscillation properties then you have

$$\lambda_{1,2} = \frac{-d \pm \sqrt{d^2 + 4ac}}{2} = -\zeta\omega_n\pm\omega_n\sqrt{\zeta^2-1} = -\frac{\omega_n}{\sqrt{2}}\pm\omega_n\sqrt{-0.5}$$ where $\zeta$ is the damping factor you want to match to 0.707 and $\omega_n$ is the natural frequency which will give you the mode shapes.

• Right, so my question is how to optimize $d$ if I don't know $a,b,c$. I'm leaning towards something like adjusting $d$ by $\frac{d(\dot x)^2}{dd}$. – Neil G Feb 3 '18 at 9:22
• @NeilG $b$ term is not relevant for the internal dynamics of $x,y$ and defines how much external forcing is affecting. For the other terms you now have a relationship between $d$ and $a,c$ that you would like to place the poles. So you can (semi)optimize such that the poles always fall inside a cone (or any region) in the complex plane. Note that this is independent from $x(t),y(t)$. – percusse Feb 3 '18 at 9:25
• It seems like it's only independent of $x(t), y(t)$ if you know $a,c$, which I don't. – Neil G Feb 3 '18 at 9:39
• @NeilG If you have no knowledge of a, c then you need to perform some sort of an system identification experiment. Otherwise there is no way to select $d$ – percusse Feb 3 '18 at 9:44
• @percusse, You haven't really said anything more than what I have already covered except the goal of $\zeta = \sqrt{1/2}$, instead of the critical value of $\zeta = 1$ which the OP is seeking. What is the rationale for that? – Cedron Dawg Feb 3 '18 at 14:14
{"extraction_info": {"found_math": true, "script_math_tex": 0, "script_math_asciimath": 0, "math_annotations": 0, "math_alttext": 0, "mathml": 0, "mathjax_tag": 0, "mathjax_inline_tex": 1, "mathjax_display_tex": 1, "mathjax_asciimath": 0, "img_math": 0, "codecogs_latex": 0, "wp_latex": 0, "mimetex.cgi": 0, "/images/math/codecogs": 0, "mathtex.cgi": 0, "katex": 0, "math-container": 0, "wp-katex-eq": 0, "align": 1, "equation": 0, "x-ck12": 0, "texerror": 0, "math_score": 0.9241424798965454, "perplexity": 301.7873277717283}, "config": {"markdown_headings": true, "markdown_code": false, "boilerplate_config": {"ratio_threshold": 0.18, "absolute_threshold": 10, "end_threshold": 15, "enable": true}, "remove_buttons": true, "remove_image_figures": true, "remove_link_clusters": true, "table_config": {"min_rows": 2, "min_cols": 3, "format": "plain"}, "remove_chinese": true, "remove_edit_buttons": true, "extract_latex": true}, "warc_path": "s3://commoncrawl/crawl-data/CC-MAIN-2021-17/segments/1618038468066.58/warc/CC-MAIN-20210418043500-20210418073500-00092.warc.gz"}