None Week03A_Modeling_MathematicalTools

Table of Contents

Challenge

This notebook's challenge is a continuation of the challenge you saw in Week 2 Friday. In this assignment, you will continue on a path towards building a suitably predictive mathematical model of your friend's hydroelectric dam system, whose behavior is represented in the simulation below.

In today's portion of this challenge, you will use data that you collected in the last assignment to build a mathematical model of the dam's response in the form of a differential equation. You will characterize this differential equation to provide more insights about how the dam acts, and imporove the predictive power of your model.

Model Development steps

We will continue to work primarily inside the first two "blocks" of the ME480 disciplined process, with diversions into controller design as needed to help us with model development:

image-2.pngUntil now, we have treated the "model development" step as if it is monolithic, but in reality, model development consists of two distinct operations:

Model Scoping

"Scoping" a mathematical model for a system means deciding what its inputs and outputs should be, and what degree of complexity the model should have. If it is a physics-based model, this complexity is often dictated by how many unique elements of the system you will model. If the model is purely mathematical, the complexity of the model may simply be dictated by the order and linearity of the differential equation(s) you choose to use to fit the model's behavior.

Model Construction

In this step, a model is actually constructed to predict the behavior of your system. For empirical, purely mathematical models, this can involve "fitting" an equation to data. For physics-based models, this can involve combining mathematical models of individual components to produce an equation or set of equations that represents the system's behavior based on known or measurable physical quantities.

In this notebook, we will take a closer look at the above two steps. Given the prerequisites for ME480, the tools presented in the sections below may look familiar. If so, that's great! However, it will still be helpful to study he following sections through the lens of the ME480 disciplined process.

Model Development: Data-Driven System Models using Differential Equations

Dynamic systems that store and dissipate energy are often well-approximated using differential equations. Sometimes, we get lucky and can use a linear differential equation to adequately model our system's dynamics. The solutions to linear differential equations satisfy the principle of superposition, and allow us to use a rich set of tools to understand our system's behavior.

When a system cannot adequately be modeled using a linear differential equation, it can often be linearized and treated like a linear system for a limited range of inputs. If it cannot be linearized, numerical integration is still possible as a way to understand how the system will behave.

The following sections focus on the development of linear models for system behavior using a data-driven approach. We will pause at times to connect these models to systems' physics, but we will leave a deep dive into physics-based modeling for a later notebook.

Model Order

In general, systems we encounter "in the wild" would require and infinite order model to capture all of the system's behavior exactly. Moreover, nearly all real systems are nonlinear. However, many times a system can be approximated by a lower-order, linear differential equation. When we say "order," we mean: how many derivatives are in the model's governing differential equation? In choosing an order for our model, an engineer's job is to count the number of significant, independent energy-storing elements in the system. That number gives us the minimum model order we can use to represent the system's dynamics.

When we look at a dataset from an experiment performed on a system, if the data "looks first order," it is often a clue that the system only has one significant independent energy storing element. If the data "looks second order," we might suspect that two of its energy storing elements are independent and significant. If the model displays a combination of first and second order behavior, it may be time to consider that the model we construct should be greater than second order. Developing an expectation for the model's order is an important step in scoping and constructing an appropriate model for a dynamic system.

In order to make this type of determination about model order, and work towards an informed model of our system, it's important to understand how first order, second order, and higher-order differential equations act.

Solving Differential Equations: the method of undetermined coefficients

The first way most students learn to solve linear, constant-coefficient differential equations is called the "method of undetermined coefficients." This is the first method we will use here, and its steps are as follows.

Given a linear differential equation of order $N$ with constant coefficients in the form:

$$a_N \frac{d^N y}{dt^N} + a_{N-1} \frac{d^{N-1} y}{dt^{N-1}}+ a_{N-2} \frac{d^{N-2} y}{dt^{N-2}} + \cdots + a_1 \frac{dy}{dt} + a_0 y = u(t)$$

Which can be written compactly as:

$$\sum_{n=0}^{n=N} a_n \frac{d^n y}{dt^n} = u(t)$$

We can say that solutions to this differential equation will satisfy the principle of superposition because the differential equation itself is linear. Therefore, we can separate the problem into finding a "homogeneous solution" where $u(t)=0$, and then adding to this a particular solution in which $u(t)$ is equal to a known function.

Because exponential functions (including complex exponential functions) are the only known functions for which their derivatives are scaled versions of themselves, we know our solution $y(t)$ must take the form of an exponential function. Therefore, we can write the "characteristic equation" for our differential equation by substituting $y(t) = e^{pt}$ into the equation's homogeneous form, where $p$ is an unknown scale factor in the exponential. When we do this, each term $\frac{d^n y}{dt^n}$ becomes $a_n p^n e^{pt}$. By canceling the common exponential term $e^{pt}$ from all terms in the equation, we are left with a polynomial (algebraic) function: $$\require{cancel}$$ $$ a_N p^N \cancel{e^{pt}} + a_{N-1}p^{N-1} \cancel{e^{pt}} + a_{N-2}p^{N-2} \cancel{e^{pt}} + \cdots + a_1 p \cancel{e^{pt}} + a_0 \cancel{e^{pt}} = 0$$ This can be written in compact form as: $$\sum_{n=0}^{n=N} a_n p^n = 0$$

The $N$ solutions to this polynomial equation are called the "characteristic roots" or "eigenvalues" of the system. Using the eigenvalues, the known input function $u(t)$ and $N$ known initial conditions, we can solve the differential equation using the following steps:

  1. Find the homogeneous solution $y_h(t)$ to the differential equation ($u(t)=0$), in terms of the system's known characteristic roots (eigenvalues) and any undetermined integration coefficients. Solve for the system's eigenvalues using the characteristic equation, but leave the integration coefficients unknown.
  2. Find the particular solution $y_p(t)$ to the differential equation, where the system is subjected to a particular input ($u(t)\neq0$)
  3. Find the total response of the system by summing the homogeneous and particular solutions.
  4. Solve for all unkown integration coefficients using the equation's $N$ known initial conditions, its $N$ characteristic roots (eigenvalues), and the known input function $u(t)$.

Stability of Linear, Constant-Coefficient Differential Equations

Because solutions to linear, constant-coefficient differential equations have solutions that include an exponential function $y(t) = e^{pt}$ as a multiplier, where an eigenvalue $p$ can be either a real number or a complex conjugate pair, we can say that a differential equation model for a system is stable (reaches a steady-state value) if and only if the real parts of all eigenvalues $p$ are strictly less than zero.

The reason for this is intuitive-- consider that if $\alpha$ is some infinitesimal positive number $e^{\alpha t}$ will have a final value of $\infty$ as $t\rightarrow \infty$, which is not a "bounded output." conversely, $e^{-\alpha t}$ will decay to zero, meaning that the models output $y$ will stop changing as $t\rightarrow \infty$.

Scoping a linear differential equation model using collected data

Scoping a differential equation model for use to describe a collected dataset means:

  1. Determining the magnitude and form of the input applied to the system when the data were collected. Was it a step response test? How large was the change in input? At what time did it occur?
  2. Collecting and plotting data of the output (and input, if possible) during the test
  3. Confirming that the system in question appears to be stable and is approximately linear for ranges of inputs you're interested in
  4. Recognizing the behavior of the system's output as closest to a first, second, or higher-order linear differential equation

Once the model is scoped, you are ready to construct it by comparing the behavior of your real system with a differential equation model.

Constructing a model using collected data

To construct a model from collected data, you will need fit a differential equation to your collected data by characterizing the collected data in terms of the differential equation model you chose. Once you do this, you can use the characteristics of the collected response to find unknown parameters in your differential equation model.

Characterizing a response in terms of a differential equation model

Characterizing (or fully describing) a response in terms of a particular differential equation model requires:

  1. Finding the approximate eigenvalues your differential equation model will need in order to match your system's transient response.
  2. Finding the steady state gain your differential equation model will need to have in order to match your system's steady-state gain.

Once these two things are known, you can equate terms in order to build a complete numerical differential equation that fits your collected data.

First Order Linear Differential Equations

A generic first-order linear differential equation with constant coefficients has the form:

\begin{equation} \sum_{n=0}^{n=1} a_n \frac{d^n y}{dt^n}=a_1\dot{y} + a_0y= u(t) \end{equation}

Because this equation is linear and satisfies the principle of superposition, the solution $y(t)$ to this differential equation can be separated into two parts:

$$y(t) = y_h(t) + y_p(t)$$

Where $y_h(t)$, or the "homogeneous solution" is the response of the system when $u(t)=0$ and $y_p(t)$, the "particular solution," is the response of the system to some specific input $u(t)$.

Characteristic equation and eigenvalue

As with all linear differential equations, the stability of a first order system is governed by its characteristic equation. Using the method of constructing the characteristic equation explaned above, we find:

\begin{equation} a_1 p + a_0 = 0 \end{equation}

Solving this algebraic equation gives us the system's eigenvalue. In german "eigen" means "own," and this value, which you will recognize as the multiplier in the exponent of the homogeneous solution, is called the system's "own" value because it is a feature of the system's solution that never changes as long as the differential equation remains the same in homogeneous form.

For a first-order system, the eigenvalue or characteristic root is:

\begin{equation} p = -\frac{a_0}{a_1} \end{equation}

The system is stable if $p<0$.

Free Response

A system's "free response" refers to the special case in which $u(t) = 0$, so no particular solution must be found. In order for the time-domain response $y(t)$ of the system to be a solution to our first order differential equation for any arbitrary nonzero initial condition $y(0) = y_0$, the homogeneous solution must have the form

\begin{equation} y_h(t) = Ce^{-\frac{a_0}{a_1}t} \end{equation}

This is because exponential functions (and trigonometric functions, which can be written as complex exponentials) are the only known functions for which the function's derivative is a scaled version of itself!!

If the initial condition $y(0) = y_0$ is known, the free response can be determined:

\begin{equation} y(0) = Ce^{-\frac{a_0}{a_1}0} = y_0 \end{equation}\begin{equation} C = y_0 \end{equation}
\begin{equation} y(t) = y_0e^{-\frac{a_0}{a_1}t} \end{equation}

Step Response

Now consider if the system has a step input of magnitude $U$, or $u(t) = Uu_s(t)$ applied, where $u_s(t)$ is the unit step function:

\begin{equation} U(t) = Uu_s(t) \end{equation}

Then the system's differential equation can be written:

\begin{equation} a_1\dot{y} + a_0y= U u_s(t) \end{equation}

The equation's homogeneous response still has the form:

$$y_h(t) = Ce^{-\frac{a_0}{a_1} t}$$

But we will leave the constant $C$ undetermined until we find the system's particular solution.

Because the homogeneous solution $y_h(t)$ decays to 0 as $t\rightarrow\infty$ (assuming $a_1 >0$ and $a_0>0$), we can find the particular solution $y_p(t)$ by considering the system's behavior at steady state. Because $u(t)$ is a constant at steady state for a step input, this tells us that the system's particular solution is also a constant! \begin{equation} y_p(t) = y_p = y_{ss} \end{equation}

Where we use $y_{ss}$ to represent the steady-state value of our equation $y(t)$. What is this steady state value? We need only to substitute back into our differential equation to find out. $\require{cancel}$

\begin{equation} \cancel{a_1\dot{y}} + a_0y= U \end{equation}

where $\dot{y}\rightarrow 0$ as $t\rightarrow \infty$. Substituting $y_{ss}$ in for $y$ as $t\rightarrow \infty$, we get:

\begin{equation} a_0y_{ss}= U \end{equation}

Solving for $y_{ss}$ yields:

\begin{equation} y_{ss} = \frac{U}{a_0} = y_p(t) \end{equation}

Now, combining the homogenous and particular solutions to obtain our total solution $y(t)$ yields:

\begin{equation} y(t) = y_p(t) + y_h(t) = \frac{U}{a_o} + Ce^{-\frac{a_0}{a_1}}t \end{equation}

Solving for $C$ knowing $y(0) = y_0$ gives us:

\begin{equation} C = y_0 - \frac{U}{a_o} \end{equation}

The complete response or a first-order linear differential equation to a step input of magnitude $U$ is: \begin{equation} y(t) = \frac{U}{a_0}(1-e^{-\frac{a_0}{a_1}t})+ y_oe^{-\frac{a_0}{a_1}t} \end{equation}

NOTE: We could have found the same result by looking at the zero-initial-condition step response and the initial-condition free response, and simply adding the solutions together! This is another important result of the principle of superposition.

Characterizing First Order Systems

Knowing the form of a first-order system's free and step responses is helpful for recognizing "approximately first order" behavior in a real system and deciding that a first order model scope may be appropriate. But in looking at data we think might be first order, we can also get a lot of information about what our mathematical model might look like by relating generic concepts like stability, settling time, and steady-state gain to mathematical features of a first order model.

Time Constant

For a first order system, the system's single eigenvalue $p=-\frac{a_0}{a_1}$. If we define the time constant to be $\tau = \frac{a_1}{a_0} = -\frac{1}{p}$, then we can use it to characterize the shape of a first order system.

When $t = \tau$, then the differential equation's homogeneous solution is \begin{equation} y_h(t) = y_0e^{-\frac{t}{\tau}} = y_oe^{-1} = 0.367879y_0 \end{equation}

This means that if $y(\tau) = 0.368 y_0$, then approximately 63.2% of the change from $t_0$ to $t_{ss}$ has occurred when $t=\tau$. For step responses, this means that the time constant can be pulled off of a plot by finding the time $\tau$ at which $y(t) = 0.632(y_{ss}-y_0)$.

image.png

Although a first order equation's homogeneous response will technically NEVER reach a value of exactly 0, and its step response will never reach a value of exactly $y_ss$, we often consider the transients "dead" and they equation "at steady state" after the settling time has passed.

Settling Time

To determine the settling time for a first-order model, we need only to refer to the definition of settling time, and then substitute this concept into your solution $y(t)$ for a first-order response. Because the exponential portion of the response is responsible for the equation's transient behavior, and because this transient behavior is due to the homogeneous response, we only need to consider the homogeneous solution when determining a system's settling time.

For the 2% settling time of a first order differential equation, this means determining when $y_h(t)<.02y_0$. Knowing that $y_h(t) = y_0e^{-\frac{t}{\tau}}$, we can make a substitution $t = t_{s,2\%}$ as we define $t_{s,2\%}$ as our "two percent settling time":

\begin{equation} 0.02y_0 = y_0e^{-\frac{t}{\tau}} \end{equation}

canceling the repeated factor $y_0$ appearing on both sides of the equation, and taking the natural logarithm of both sides, we see that:

\begin{equation} ln(.02) = ln(e^{-\frac{t_{s,2\%}}{\tau}}) = -3.91 \end{equation}

This means that:

\begin{equation} t_{s,2\%} = 3.91\tau \end{equation}

This is usually rounded up, and we typically write:

\begin{equation} t_{s,2\%} \approx 4\tau \end{equation}

Steady State Gain

To determine steady-state gain for a first-order model's step response, we simply need to look to the definition of steady state gain and apply it to our equation. Because we are talking about steady state behavior, we only need to consider the particular solution to the differential equation when the system is subjected to a step response of magnitude $U$:

$$a_0 y_p = U$$

as we learned above, this means that the steady state value of $y$ can be written:

$$y_p = y_{ss} = \frac{U}{a_0}$$

The steady state gain is defined as the ratio of the change in output to the change in input for a system. Because our model is a linear differential equation, we know that superposition applies, so it's possible that the total input $U$ at steady state is the result of more than one step input. In fact, if the system was moving at the beginning of a test, it must have some prior input that caused it to move! Therefore, considering the possibility of a "prior" step input that created an initial condition $y_0$ in our equation, we typically write:

$$\Delta y_{ss} = y_{ss}-y_0 = \frac{U-U_0}{a_0}$$

This allows us to write an equation for steady state gain using the first-order differential equation's $a_0$ term.

$$K_{ss} = \frac{y_{ss}-y_0}{U-U_0} = \frac{\Delta y_{ss}}{\Delta u_{ss}} = \frac{1}{a_0} $$

This concludes our study of first-order models for now-- in the sections below, second order models are explored, but there is slightly less detail in the derivation steps. The process for obtaining the results below is the same mathematical process as we used for first-order models. More details about the derivation of any piece of the material below can be found in a differential equations textbook or by asking your instructor.

Second-Order Linear Differential Equations

A linear, constant-coefficient second-order differential equation has the form:

$$\sum_{n=0}^{n=2} a_n \frac{d^n y}{dt^n}=a_2 \ddot{y} + a_1\dot{y} + a_0 y = u(t)$$

Because this equation is linear and satisfies the principle of superposition, the solution $y(t)$ to this differential equation can be separated into two parts:

$$y(t) = y_h(t) + y_p(t)$$

Where $y_h(t)$, or the "homogeneous solution" is the response of the system when $u(t)=0$ and $y_p(t)$, the "particular solution," is the response of the system to some specific input $u(t)$. In order to fully determine a second order system's response, it is vital to know the initial conditions $y(0)$ and $\dot{y}(0)$ as well as the system's eigenvalues.

Eigenvalues and characteristic equation

As with a first order system, the characteristic equation for a second order system is obtained by substituting $\frac{d^n y}{dy^n}$ with the dummy variable $p^n$ into the equation's homogeneous form with $u(t)=0$. For a second order system, this yields:

$$a_2 p^2 + a_1 p + a_0 = 0$$

This equation can be solved using the quadratic formula to obtain the system's eigenvalues, also called the system's poles or characteristic roots. The quadratic formula will always yield two solutions to the characteristic equation for a second-order system.

$$p_1,p_2 = \frac{-a_1 \pm \sqrt{a_1^2 - 4 a_2 a_0} }{2a_2}$$

This equation leaves three possibilities:

  • The roots $p_1,p_2$ are the same number. This system is referred to as having "repeated roots."
  • The roots $p_1,p_2$ are two different real numbers
  • The roots $p_1,p_2$ are complex conjugates

In all cases, stability of the system depends on the real parts of $p_1,p_2$ both being strictly negative. For an unstable system, the second-order equation will result in an infinite value at steady state. For a stable system, transient behavior depends on whether the eigenvalues are real and repeated, real and distinct, or complex conjugates.

Free Response: general form for real, distinct characteristic roots

Like a first-order differential equation, the solution to a generic second order differential equation's homogeneous response will take the form of an exponential, but will contain two exponential terms rather than one, since there are two derivatives in the equation:

$$ y_h(t) = C_1 e^{p_1 t} + C_2 e^{p_2 t}$$

for the case of real, distinct roots $p_1,p_2$, we can use the known initial conditions $y(0)= y_0,\dot{y}(0) = \dot{y}_0$ and substitute them into the differential equation to find that:

$$y_h(t) = \left(y_0 - \frac{\dot{y}_0 - p_1 y_0}{p_2 - p_1}\right) e^{p_1 t} + \left(\frac{\dot{y}_0 - p_1 y_0}{p_2 - p_1}\right)e^{p_2 t} $$

Free Response: general form for complex conjugate roots

For complex conjugate roots, we can still write the homogeneous solution in the same form as for real, distinct roots:

$$ y_h(t) = C_1 e^{p_1 t} + C_2 e^{p_2 t}$$

However, we can now also write $p_1 = \sigma - j\omega_d$ and $p_2 = \sigma+j\omega_d$, since the eigenvalues are complex conjugates. By using Euler's Formula, $e^{-j \omega_d} = cos\omega_d + j \sin \omega_d$, and recognizing that the constants $C_1$ and $C_2$ must also be complex conjugates for the solution $y_h(t)$ to be real for all times $t$, we can obtain an equation the total free response for the system in terms of the eigenvalues' real and imaginary parts.

$$y_h(t) = e^{\sigma t}\left( C_3cos(\omega_d t) + C_4sin(\omega_d t) \right)$$

Here, $C_3$ and $C_4$ are real constants, and we can find them using the known initial conditions $y(0) = y_0$ and $\dot{y}(0) = \dot{y}_0$. Every term in the above equation is real, because we know for a fact that the response of our differential equation will be real, so all terms multiplied by the complex number $j$ as a result of the Euler formula expansion must cancel.

Solving for the free response means solving for the coefficients $C_3$ and $C_4$ in the homogeneous equation using the initial conditions and the system eigenvalues $p_1,p_2 = \sigma \pm j \omega_d$.

$$y_h(t) = e^{\sigma t}\left( y_0 cos(\omega_d t) + \frac{\dot{y}_0 - \sigma y_0}{\omega_d}sin(\omega_d t) \right)$$

Free Response: general form for repeated roots

We can think of the case of repeated roots, where $p_1=p_2$, as a special case of the complex conjugate roots case where $p_1,p_2 = \sigma \pm j\omega_d$, where we shrink the imaginary part $\omega_d$ of the complex conjugate pair to 0 to get two roots $p_1,p2 = \sigma$. Using the form of the free response from the complex conjugate roots case, we see that:

$$y_h(t) = e^{\sigma t}\left( C_3cos(\omega_d t) + C_4sin(\omega_d t) \right) =e^{\sigma t}\left( y_0 cos(\omega_d t) + \frac{\dot{y}_0 - \sigma y_0}{\omega_d}sin(\omega_d t) \right)$$

In the limit as $\omega_d \rightarrow 0$, we get a division by zero in the $C_4$ term! Let's look at that more closely:

$$\lim_{\omega_d \rightarrow 0} \frac{\dot{y}_0 - \sigma y_0}{\omega_d} = \frac{\dot{y}_0 - \sigma y_0}{0}$$

Not good!! To figure out what to do next, let's try the limit of the entire $C_4 \sin(\omega_d t)$ term.

$$\lim_{\omega_d \rightarrow 0} \frac{\dot{y}_0 - \sigma y_0}{\omega_d} \sin(\omega_d t) = \frac{\dot{y}_0 - \sigma y_0}{0} \cdot 0$$

Still a problem! But let's try using L'Hopital's rule for limits. Pulling the known constant $\dot{y}_0 - \sigma y_0$ out of the limit, we find:

$$(\dot{y}_0 - \sigma y_0)\left[\lim_{\omega_d \rightarrow 0} \frac{\sin(\omega_d t)}{\omega_d } = \lim_{\omega_d \rightarrow 0} \frac{\frac{d}{d\omega_d} \sin (\omega_d t)}{\frac{d}{d\omega_d}\omega_d} =\lim_{\omega_d \rightarrow 0} \frac{t \cos(\omega_d t)}{1} = t\right]$$

Therefore, the free response of a second order linear differential equation when the equation has two repeated eigenvalues $p_1,p_2 = \sigma $ is:

$$y_h(t) = e^{\sigma t}\left( y_0 cos(\omega_d t) + (\dot{y}_0 - \sigma y_0)t \right)$$

Substituting in $\omega_d = 0$ yields a further simplification, since $\cos(0) = 1$.

You may have seen this in a differential equations course in a general form that looks something like this:

$$y_h(t) = e^{\sigma t}(C_5 + t C_6)$$

Where $C_5$ and $C_6 $ can be found using the initial conditions and the eigenvalue $\sigma$

$$y_h(t) = e^{\sigma t}(y_0 + t(\dot{y}_0 - \sigma y_0))$$

Step Response with Zero Initial Conditions: Real, Distinct Roots

Limiting our discussion to step responses with $y(0) = 0, \dot{y}(0)=0, u(t) = U u_s(t)$, we can find the particular solution by considering the equation's behavior at steady state:

$$y_p(t) = y_{ss} = \frac{U}{a_0}$$

Combining the particular and homogeneous solutions for real distinct roots yields:

$$ y(t) = C_1 e^{p_1 t} + C_2 e^{p_2 t} + \frac{U}{a_0}$$

With zero initial conditions, the constants $C_1$ and $C_2$ can be found as:

$$ C_1 = \frac{U}{a_0}\left( 1 - \frac{p_1}{p_2-p_1}\right)$$$$ C_2 = \frac{p_1}{p_2 - p_1}\frac{U}{a_0}$$

Which yields the following equation for the zero-initial-condition step response:

$$ y(t) = \frac{U}{a_0}\left( \frac{p_2}{p_1-p_2} e^{p_1 t} + \frac{p_1}{p_2 - p_1} e^{p_2 t} + 1\right)$$

Step Response with Zero Initial Conditions: Complex Conjugate Roots

Limiting our discussion to step responses with $y(0) = 0, \dot{y}(0)=0, u(t) = U u_s(t)$, we can find the particular solution by considering the equation's behavior at steady state:

$$y_p(t) = y_{ss} = \frac{U}{a_0}$$

Combining the particular and homogeneous solutions for real distinct roots yields: $$y(t) = e^{\sigma t}\left( C_3cos(\omega_d t) + C_4sin(\omega_d t) \right) + \frac{U}{a_0}$$

With zero initial conditions, the constants $C_3$ and $C_4$ can be found as:

$$C_3 = -\frac{U}{a_0}$$$$C_4 = \frac{\sigma U}{\omega_d a_0}$$

Which leads to the following equation for the zero-initial-condition step response:

$$y(t) = \frac{U}{a_0} \left[ 1 - e^{\sigma t}\left(cos(\omega_d t) - \frac{\sigma}{\omega_d } sin(\omega_d t) \right) \right]$$

Step Response with Zero Initial Conditions: Repeated Roots

Limiting our discussion to step responses with $y(0) = 0, \dot{y}(0)=0, u(t) = U u_s(t)$, we can find the particular solution by considering the equation's behavior at steady state:

$$y_p(t) = y_{ss} = \frac{U}{a_0}$$

Combining the particular and homogeneous solutions for real distinct roots yields:

$$y(t) = e^{\sigma t}(C_5 + t C_6) + \frac{U}{a_0}$$

With zero initial conditions, the constants $C_5$ and $C_6$ can be found as:

$$C_5 = -\frac{U}{a_0}$$$$C_6 = -\frac{\sigma U}{a_0}$$

Which yields the following equation for the zero-initial-condition step response:

$$y(t) = \frac{U}{a_0}\left[ 1 - e^{\sigma t} \left( 1 + \sigma t \right) \right]$$

Step Responses with Nonzero Initial Conditions

To find step responses of second order linear differential equations that have nonzero initial conditions can be achieved by either following the method of undetermined coefficients as above without the simplification that $y(0) = \dot{y}(0) = 0$, or by using the principle of superposition. By the principle of superposition, one could add the "free response" for a second order system to the "zero initial condition" step response to obtain the total response of the system.

This concludes our discussion of free and step responses of second order linear differential equations. As with first order systems, the behavior of second order equations can be characterized to help us recognize second order behavior in data.

Characterizing Second-Order Responses

Knowing the form of a second-order system's free and step responses is helpful for recognizing "approximately second order" behavior in a real system and deciding that a second order model scope may be appropriate. But in looking at data we think might be second order, we can also get a lot of information about what our mathematical model might look like by relating generic concepts like stability, settling time, and steady-state gain to mathematical features of a second order model. In addition to these general concepts, second-order systems with complex conjugate eigenvalues oscillate and decay if they are stable. Looking at how these systems oscillate and decay can help us infer physical insights from data that look second order.

image.png

Characteristic Equation: Forms for Overdamped, Underdamped, and Critically Damped Equations

Consider a second-order linear differential equation in "standard form": $$a_2 \ddot{y} + a_1\dot{y} + a_0 y = u(t)$$

Its characteristic equation is:

$$a_2 p^2 + a_1 p + a_0 = 0$$

After solving for the system's characteristic roots, we will know whether it is:

  1. Underdamped, with a complex conjugate pair for the eigenvalues $p_1,p_2$
  2. Critically damped, with two repeated real eigenvalues $p_1 = p_2$
  3. Overdamped, with two distinct real eigenvalues $p_1,p_2$

Each of these types of systems lends itself to a different "standard form" for the characteristic equation.

For the system to be stable, the real parts of both of its eigenvalues must be strictly less than 0.

Underdamped systems

For an underdamped system, the standard form of the characteristic equation is:

$$p^2 + 2\zeta \omega_n p + \omega_n^2 = 0$$

By equating coefficients, we can see that $2\zeta \omega_n = \frac{a_1}{a_2}$ in our original differential equation, and $\omega_n^2 = \frac{a_0}{a_2}$.

Critically Dampled systems

For a critically damped system, the characteristic equation is often factored:

$$\left(p+\frac{1}{\tau}\right)^2 = 0$$

Where $\tau$ is the "effective time constant" of the system. By equating coefficients, we can see that $\frac{2}{\tau} = \frac{a_1}{a_2}$ in our original equation, and $\frac{1}{\tau^2} = \frac{a_0}{a_2}$ in our original equation.

Overdamped systems

For an overdamped system, the characteristic equation is often factored:

$$\left(p+\frac{1}{\tau_1}\right) \left(p+\frac{1}{\tau_2}\right) = 0$$

The "effective time constants" $\tau_1,\tau_2$ can be found by equating coefficients with our original equation as $\frac{1}{\tau_1} + \frac{1}{\tau_2} = \frac{a_1}{a_2}$ and $\frac{1}{\tau_1 \tau_2} = \frac{a_0}{a_2}$.

Unlike first-order systems, we cannot, in general, say that the response of an overdamped or critically damped second-order system is 63.2% finished with its transient behavior at one time constant $\tau$, but for overdamped systems, this approximation becomes better and better as the eigenvalues $p_1$ and $p_2$ become more and more different, with one decaying "much more quickly" than the other.

Damping Ratio, Natural Frequency, and Damped Natural Frequency for underdamped systems

Looking at the standard form of the characteristic equation for an underdamped second-order equation:

$$p^2 + 2\zeta \omega_n p + \omega_n^2 = 0$$

We can see that the quadratic formula yields solutions $p_1,p_2$ of:

$$p_1,p_2 = - \zeta \omega_n \pm \omega_d j $$

Where $\omega_d = \omega_n \sqrt{1-\zeta^2}$ is called the "damped natural frequency" of the system. $\omega_d$ is the same term we used in the section above on homogeneous and step responses for underdamped systems, and it represents the imaginary part of the eigenvalue. It tells us the frequency in radians per second at which the equation oscillates. The "natural frequency" $\omega_n$ tells us how fast the equation would oscillate in the absence of damping, or with a damping ratio of $\zeta = 0$. Note that when $\zeta=0$, the real part of the complex conjugate pair disappears, and the system oscillates at $\omega_d = \omega_n$.

To see how damping ratio, natural frequency, and damped natural frequency are related, see the example code below!

In [2]:
%let's plot the step response for a prototype second order system with a few damping ratios.

zeta = .2:.2:2;%damping RATIO
wn = 10*2*pi;%natural frequency
%get us a time vector!
t = linspace(0,.5,1000);
%initialize our solution matrix
y = zeros(length(zeta),1000);
%loop through zetas applying the correct solution.
%note that the form of the solution here is different but equivalent to the forms derived above (CHECK THEM IF YOU LIKE!)
for k = 1:length(zeta)
    if zeta(k)<1
        y(k,:) = 1-1/sqrt(1-zeta(k)^2).*exp(-zeta(k)*wn*t).*sin(sqrt(1-zeta(k)^2).*wn*t+acos(zeta(k)));
        legendInfo{k} = ['\zeta = ' num2str(zeta(k))];
    elseif zeta(k)==1
        y(k,:) = (1-exp(-zeta(k)*wn.*t)-zeta(k)*wn.*t.*exp(-zeta(k)*wn.*t));
        legendInfo{k} = ['\zeta = ' num2str(zeta(k))];
    else
        a = abs(-zeta(k)*wn + wn*sqrt(zeta(k)^2-1));
        b = abs(-zeta(k)*wn - wn*sqrt(zeta(k)^2-1));
        y(k,:) = a*b*1/(a*b)*(1-b*exp(-a*t)/(b-a)+a*exp(-b*t)/(b-a));
        legendInfo{k} = ['\zeta = ' num2str(zeta(k))];
    end
end

figure
plot(t,y)
legend(legendInfo)
xlabel('Time (s)')
ylabel('y(t)')

Finding an underdamped system's damped natural frequency and damping ratio from a plot

To find $\omega_d$ and $\zeta$ from a plot, which will allow you to find both of the system's eigenvalues, first note that the damped natural frequency can be found using the time $T$, or period, between peaks in an oscillatory response:

$$\omega_d = \frac{2\pi}{T}$$

The system's damping ratio $\zeta$ can be found using the log decrement formula.

\begin{equation} \zeta = \frac{\frac{1}{n-1}\left(ln\frac{y_1}{y_n}\right)}{\sqrt{4\pi^2+\left(\frac{1}{n-1}\left(ln\frac{y_1}{y_n}\right)\right)^2}} \end{equation}

Where the definitions of $y_1,y_2,\ldots,y_n$ are given by the following figure: image.png

Settling Time

For lightly damped underdamped systems (with small damping ratio $\zeta$, the settling time is approximated similarly to the method used for first-order systems. Because the eigenvalues' real parts are located at $Re(p) = -\zeta \omega_n = \sigma$, the step response of a second order system is bounded by an exponential function $e^{\sigma t} = e^{-\frac{t}{\tau_{eff}}}$, where $\zeta \omega_n = \frac{1}{\tau_{eff}}$ can be thought of as an "effective time constant" for the system. This means that the 2% settling time can be approximated as:

$$t_{s,2\%} \approx \frac{4}{\zeta \omega_n} = 4 \tau_{eff}$$

For overdamped systems with large separations between $p_1$ and $p_2$, the slower (larger) effective time constant can be used to compute settling time. For critically-damped or nearly critically-damped systems, the settling time can be computed analytically by finding the point at which the solution never leaves the 2% bounds around its total change $\Delta y_{ss} = y_{ss}- y_0$.

Steady State Gain

Still the ratio of change in input over change in output, the steady state gain for any standard-form second-order system in response to a step input can be computed as:

$$K_{ss} = \frac{y_{ss} - y_0}{U_{ss} - U_0} = \frac{\Delta y_{ss}}{\Delta u_{ss}} = \frac{1}{a_0}$$

While these aren't the only ways to characterize a second order system, these are the most common tools we'll use in ME480.

Higher-Order Linear Differential Equations and the method of undetermined coefficients

First and second order systems are nice when they describe our real system's behavior. They are also nice because finding the response of a higher-order (greater than 2nd order) system can be accomplished using the method of undetermined coefficients-- higher-order systems can only display combinations of first and second-order behavior, because the characteristic equation can only have either real or complex conjugate eigenvalues!

This doesn't change the fact that using the method of undetermined coefficients to solve 3rd and higher order systems can be tedious. We will develop more sophisticated, more efficient tools to deal with higher-order systems as we need them.

Exercise 1: Due at Midnight on Wednesday Sept. 15

Model Development: Develop a first or second order linear model of the Dam's behavior

Using your code from Week 2 Friday, collect a dataset with the dam overflowing, the turbine at rest, and the penstock valve stepping from closed to 25% open, effectively stepping the input pressure at the penstock valve from zero to some nonzero value (because the valve begins closed, and then suddenly opens to let water through). Collect data for 5 seconds. Upload your saved data to this assignment's folder so that you can load it into Octave for plotting.

Scope and then construct a first or second order model (as appropriate) to your dataset that describes the relationship between the dam's hydrostatic pressure and the voltage at the house.

Present any hand calculations that support your model development and/or fitting, along with any MATLAB code that supports your efforts. Include your model's differential equation in your final answer, along with the differential equation's time-domain solution for $V(t)$ for a step change in dam height.

Your MATLAB code must show a properly formatted and labeled plot that shows the fit between your model and the data.

YOUR ANSWER HERE

In [3]:
% YOUR CODE HERE
error('No Answer Given!')

error: No Answer Given!

Model Validation: How well does the model predict the dam's behavior?

Validate your model by using it to predict house voltage in response to a different input hydrostatic pressure. Accomplish this by (manually) draining the reservoir to approximately half full, and collecting data with the reservoir beginning at th new lower height, with the turbine stopped and the valve opening from closed to 25%. In the Markdown cell below, explain whether you think your model appropriately predicts the dam's behavior under this new set of conditions and why. Be specific by using the concepts from this notebook.

YOUR ANSWER HERE

In [4]:
% YOUR CODE HERE
error('No Answer Given!')

error: No Answer Given!