None W04A_Modeling_Nonlinear_Numerical

Table of Contents

Challenge

This notebook's challenge is a continuation of the challenge you saw in Week 3 B. 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. All of the available variables and operation are the same as they were in Wednesday's assignment.

Simulation Data Format

The data file's structure is the same as it was for your last assignment::

column 1 2 3 4 5 6
data time gauge height (ft) Pressure at penstock valve (kPa) turbine speed (RPM) generator volts house volts

As you may recall from Week02B, one of the original goals of this challenge was to develop a system to regulate the voltage at the house to make sure that the charging station has the proper voltage.

As we continue towards that larger goal, today's challenge involves investigating a possible strategy for controlling the turbine's speed. Because the voltage at the house appears directly related to turbine speed, this seems like a good next step. We will leave the house disconnected from the turbine for today's challenge, and tackle the problem of controlling turbine speed first without the added complexity of the electrical system.

The type of control we will focus on in ME480 is called classical feedback control. In classical feedback control, information about a system's output is "fed back" to the controller, which may be an Arduino or similar microprocessor, in order for it to make "informed decisions" about what to do next. A graphical representation of this paradigm is shown as a block diagram below.

image.png

In the diagram above, the controller implements a control law to turn the error between $\Omega_d$, or the "desired angular speed of the turbine," into action that it takes in order to influence the turbine. If the system could do this "on its own," it could automatically account for variations in the dam height due to rainfall or drought, keeping the turbine spinning at the correct speed all the time!

The most obvious action the dam system could take to influence the turbine's speed is to change the penstock valve position to control the fluid resistance of the valve. Closing the valve slightly could increase the resistance, slowing the turbine down by providing it less fluid power. Opening the valve slightly would decrease the valve's resistance, speeding the turbine up.

By making "decisions" about how much to open or close the valve based on the error $e$ between desired and actual turbine speed, it is possible for the controller to regulate the turbine's speed automatically. This implies that the arrow representing the input to the dam system should be valve position $u_v$, as shown below.

image-3.png

Looking at this diagram, one can see that the overall input to the system is desired turbine speed, and the output of the overall system is actual turbine speed. With a good design for the controller, these two should match fairly well, no matter what happens. Looking at the individual "pieces" of this diagram, we can see that the "controller" should have an input of error and an output of valve position, and the "dam system" should have an input of valve position and an output of turbine speed.

But wait... in your last assignment, for which you developed a physics-based model, $u_v$ was not the input! $P_{in}$, or the hydrostatic pressure from the dam was. You may have used a lumped-element physical model consistent with the following linear graph:

image-2.png

Where $P_{in}$ is the hydrostatic pressure from the dam, $R$ is the penstock and valve fluid resistance, $D$ represents an idealized turbine of displacement $D$, $b$ represents viscous friction in the turbine's shaft, and $J$ represents the turbine shaft's rotational mass moment of inertia. You may have used this linear graph to obtain a differential equation model that looked like:

$$\dot{\Omega}_{1g} = \frac{1}{J}\left[ D P_{in} - (b + D^2 R) \Omega_{1g}\right]$$

If so, you might notice that the constant $P_{in}$ (after the valve was opened) was considered our input, and the differential equation model was linear since doubling the input pressure $P_{in}$ would have doubled the output turbine speed $\Omega_{1g}$.

But now, we wish to use the valve position $u_v$ as our control "action," whereby the controller changing the valve percentage between values of $0$ and $100$ can be used to change the resistance $R$.

Before we can understand how we might design a controller that would accomplish our goals, we have to re-think our model for the dam system, and develop a physics-based model from valve position to turbine speed.

Let's assume that using some data from the valve's manufacturer, you are aware that:

$$R \approx R_{nom}(1 - .01u_v)$$

Where $R_{nom}$ is the resistance of the valve when it is just barely cracked above the 0% valve position.

If we substitute this relationship into our differential equation model above, we can see that:

$$\dot{\Omega}_{1g} = \frac{1}{J}\left[ D P_{in} - (b + D^2 R_{nom}(1-.01u_v)) \Omega_{1g}\right]$$

Breaking out terms yields:

$$\dot{\Omega}_{1g} = \frac{1}{J}\left[ D P_{in} - (b + D^2 R_{nom} - .01 D^2 R_{nom}u_v) \Omega_{1g}\right]$$

Uh-oh. If we are saying that the overall input to our system is $u_v$, and the overall output is $\Omega_{1g}$, then we can see immediately that the solution to this differential equation will not satisfy the principle of superposition, and it is therefore not linear!

THIS IS PROBABLY WHAT YOU EXPECTED Based on your initial investigation of the relationship between valve position and house voltage from the Week02B assignment. There, you found based on observation that doubling the valve position with a particular dam input pressure did not double the system output.

But how could I recognize the nonlinearity in the model, before a test was performed? A few things to look for are in the cell below.

Recognizing and Handling Nonlinearities in System Differential Equations

Differential equations are linear if they satisfy the principle of superposition. All linear, constant-coefficient differential equations can be written in the compact form:

$$\sum_{n=0}^{n=N} a_n \frac{d^n y}{dt^n} = \sum_{m=0}^{m=M} b_m \frac{d^m u}{dt^m} $$

Where $M<=N$, $u(t)$ is the equation's input, and $y(t)$ is the equation's output. We call this an "input-output" differential equation. This is slightly different than the compact form we saw here but this one allows us to consider a known input whose derivatives might also be important for describing the system's behavior. Theoretically, if the input function is known, its derivatives are also known, and the "total input" to the system could be redefined so that no input derivatives exist. However, we will consider the more general case here, and use the form above.

For this discussion, we will assume that you already have your differential equation model in input-output form.

We can recognize nonlinearities in an input-output differential equation by recognizing any deviation from the form above. Some common nonlinearities include:

  • A constant term (one that is not multiplied by an output, output derivative, input, or input derivative).
  • A term that is a nonlinear function of the input, a derivative of the output, or a derivative of the output (e.g. $\dot{y}^2$ or $\left| y \right|$).
  • Single term that includes both the input and the output (e.g. $u\cdot y$)

To handle nonlinearity in your model, you have three basic choices:

  1. Solve (or simulate) the nonlinear differential equation.
  2. Replace any nonlinear physical elements with approximately equivalent linear ones.
  3. Linearize the system about a normal operating point using the Taylor Series Expansion.

Some types of nonlinearities are very easy to handle. For example, "extra" constant terms can sometimes be treated as separate inputs, and then the "total" solution to the differential equation is the sum of the particular solutions to each. By re-imagining these terms as inputs, the system can be treated as a linear system by using superposition to find several particular solutions, and adding them together.

On the other hand, some types of nonlinearities cannot be handled using linear differential equations techniques at all. For example, the term $\left|y\right|$ has a discontinuity at $y=0$. The only option is to "solve the nonlinear differential equation (option 1)." In ME480, our only option for this type of situation is to arrive at a solution using numerical integration, which is covered below.

With any luck, however your system operates in response to a small range of inputs that are all in the neighborhood of some average known input, is stable for the range of inputs your system will experience, and thus the system's output remains in a small neighborhood of some average known output. We call the combination of an average known input and an average known output the Normal Operating Point of your system.

With a little more luck, any nonlinear terms in your differential equation can be differentiated at this normal operating point. If this is the case, you can use the Taylor Series Expansion to linearize your input-output differential equation according to the procedure below.

Linearizing an Input-Output Differential Equation using the Taylor Series Expansion

Sometimes, it is possible to linearize a general nonlinear function (like a differential equation) about a specific normal operating point (NOP), or the point where the system's output rests at steady state given a "normal" or expected steady-state input. "Linearize" in this case means turning a nonlinear differential equation into a linear one. Not all nonlinear systems can be linearized, but many can if they contain terms that are all differentiable.

Usually, this is done by approximating the nonlinear terms in a differential equation using the first two terms in the term's Taylor Series expansion. Given a generic nonlinear function $f_{NL}(y)$ and a replacement of the dependent variable $y=\overline{y}+\Delta y$, where $\Delta y$ is a small perturbation of the output variable $y$ about the system's "normal" output value $\overline{y}$, we can write the function as its infinite Taylor series:

\begin{equation} f(y) = f(\overline{y})+\left.\frac{\partial f}{\partial y}\right|_{y=\overline{y}} \left(\Delta y\right) + \frac{1}{2!}\left.\frac{\partial f}{\partial y}\right|_{y=\overline{y}} \left(\Delta y \right)^2 + \ldots \end{equation}

Which is not all that helpful except that if we just use the first term in the Taylor series, we necessarily get a linear function! This means that:

\begin{equation} f(y) \approx f(\overline{y})+\left.\frac{\partial f}{\partial y}\right|_{y=\overline{y}} \left(\Delta y\right) \end{equation}

Linearizing a generic nonlinear function of multiple variables $f_{NL} = f_{NL}(y,u)$ could similarly be approximated by the following first-order terms in the Taylor series expansion using partial fractions:

$$f_{NL}(y,u) = f_{NL}(\overline{y},\overline{u}) +\left.\frac{\partial f_{NL}}{\partial y}\right|_{y=\overline{y},u=\overline{u}} \left(\Delta y\right) +\left.\frac{\partial f_{NL}}{\partial u}\right|_{y=\overline{y},u=\overline{u}} \left(\Delta u\right)$$

In this case, we are only using the first derivative term in the Taylor Series expansion for each dependent variable. The derivatives in the expansion are constants because they are evaluated at a single point. This allows us employ the first term of the Taylor Series as a "good enough" approximation of a nonlinear function by a linear one in a small neighborhood of $y=\overline{y}, u = \overline{u}$.

For a nonlinear term $f_{NL}$ in an input-output differential equation, which might be a function of multiple variables that are not constants, we might see something more like:

$$f_{NL} = f_{NL} \left( y, \frac{dy}{dt}, \frac{d^2y}{dt^2}, \cdots , \frac{d^Ny}{dt^N}, u, \frac{du}{dt}, \frac{d^2u}{dt^2}, \cdots , \frac{d^Mu}{dt^M} \right)$$

In other words, a nonlinear term might include dependency on both the system's output, the system's input, and their respective derivatives.

If we make the assumption that the system is BIBO stable, and that the variables making up the NOP, $\overline{y},\overline{u}$ are constants, all derivatives of $\overline{y},\overline{u} = 0$. Using more first-order Taylor series terms, we can add in effects from a term's dependency on an input or output derivative by adding terms such as:

$$ + \left. \frac{\partial f_{NL}}{\partial \dot{y}} \right|_{y = \overline{y}, \dot{\overline{y}}=0,u=\overline{u}}$$

to our expansion for that term.

But all of this only addresses a single nonlinear term in our differential equation. What we need is a disciplined process to linearize an input-output model. This is given below.

Disciplined Process for Linearizing an input-output Differential Equation

We will follow a four-step process for linearizing a nonlinear input-output model, assuming you already have your system's input-output differential equation, and it is known to be nonlinear. The steps are:

  1. Determine the system's Normal Operating point, defined by a constant $u = \overline{u}$, $y = \overline{y}$. One of these two quantities must be known (either the average input or average output), and the other can be found by setting all derivatives in the differential equation to zero.
  2. Using the Taylor Series expansion, linearize all non-constant, nonlinear terms in your system's differential equation.
  3. Substitute the linearized terms back into your differential equation, substituting all remaining $y = \overline{y} + \Delta y$ and $u = \overline{u} + \Delta u$.
  4. Confirm that all constant terms now cancel out of your differential equation, leaving you with a differential equation that is linear from the new incremental input $\Delta u$ to the new incremental output $\Delta y$. This differential equation now represents your model for how the system moves in a small neigborhood around the NOP.

Example

Consider the following input-output differential equation model which was derived for some arbitrary system using a physics-based, lumped-parameter approach, where $u$ is the system's input, and $y$ is the system's output.

$$ a^2 \ddot{y} + \dot{y} e^{-\dot{y}} + uy = K$$

where $a$ and $K$ are constants.

This differential equation model is nonlinear in the following ways:

  • The nonlinear function of $\dot{y}$, which we can call $f_1(\dot{y}) = \dot{y} e^{-\dot{y}}$
  • The constant $K$ is neither an input nor an output
  • the input and output are multiplied together, in a term which we can call $f_2 = uy$

Let's follow the four-step process above to linearize this input-output differential equation model.

Step 1: Find the NOP

In order to find the system's NOP, we must assume that we either know its average input $u = \overline{u}$ or its average output $y = \overline{y}$. We will also need to assume that system doesn't stray far from these values.

We begin by setting all derivatives in the differential equation model to 0, since if we're evaluating the system at the NOP, we are assuming the values of $u$ and $y$ are constants.

$$ a^2\cdot 0 + o \cdot e^{0} + \overline u \overline{y} = K$$

This leads to the following relationship between $\overline{u}$ and $\overline{y}$:

$$\overline{y} = \frac{K}{\overline{u}}$$

Once again, we would need to either know the value of $\overline{y}$ or the value of $\overline{u}$ to calculate the total NOP.

Step 2: Use the Taylor Series to linearize each nonlinear term

Using the Taylor Series expansion's first term for each nonlinear portion of the differential equation starting with $f_1$, *evaluating each relevant term at the NOP with $y= \overline{y},u=\overline{u},\dot{y} = 0,\ddot{y}=0$, we find:

$$\require{cancel}$$$$f_1 = \dot{y} e^{-\dot{y}} \approx \left. \dot{y} e^{-\dot{y}}\right|_{\dot{y}=0} + \left. \frac{\partial}{\partial \dot{y}} \dot{y} e^{-\dot{y}} \right|_{\dot{y}=0} \Delta \dot{y}= 0 + \left. e^{-\dot{y}}\right|_{\dot{y} = 0}\Delta \dot{y} + \cancel{\left. \dot{y}^2 e^{-\dot{y}} \right|_{\dot{y}=0}\Delta \dot{y}}$$

Simplifying yields: $$f_1 \approx \Delta \dot{y}$$

For the term $f_2 = uy$, we find:

$$f_2 = uy \approx \left. uy \right|_{u=\overline{u},y=\overline{y}} + \left. \frac{\partial}{\partial y} uy \right|_{u=\overline{u},y=\overline{y}} \Delta y + \left. \frac{\partial}{\partial u} uy \right|_{u=\overline{u},y=\overline{y}} \Delta u = \overline{u}\overline{y} + \overline{u}\Delta y + \overline y \Delta u$$

Simplifying, we find that:

$$f_2 \approx \overline{u}\overline{y} + \overline{u}\Delta y + \overline y \Delta u$$

Step 3: Substitute Linearized terms back into the differential equation

Substituting $f_1$ and $f_2$ in their linearized, approximate forms back into our differential equation, and replacing any remaining $y$ and $u$ terms with $y = \overline{y} + \Delta y$, $u = \overline{u} + \Delta u$, we find:

$$ a^2 \frac{d^2}{dt^2}(\overline{y} + \Delta y) + \Delta \dot{y} + \overline{u}\overline{y} + \overline{u}\Delta y + \overline y \Delta u = K$$

We know that $\frac{d^2}{dt^2}\overline{y} = 0$ because $\overline{y}$ is a constant, so this reduces to:

$$ a^2 \ddot{ \Delta y} + \Delta \dot{y} + \overline{u}\overline{y} + \overline{u}\Delta y + \overline y \Delta u = K$$

Step 4: Confirm that all constant terms cancel, and that the differential equation is now linear in its incremental input and output

This is now a differential equation with $\Delta y$ as its output, and $\Delta u$ as its input. These values represent the small perturbations about the NOP. Using the relationship we developed in Step 1, we can see that $\overline{y}\overline{u} = K$. This means that we end up with the final differential equation:

$$a^2 \Delta \ddot{y} + \Delta \dot{y} + \overline{u}\Delta y + \overline y \Delta u = 0$$

This equation no longer has any nonlinear terms-- it only contains constant multiples of inputs, outputs, and their derivatives! We have now linearized this input-output differential equation.

Model Validation for Linearized Models

In order to determine whether your linearized model is an appropriate approximation of the nonlinear model, you can do one of two things:

  1. Compare the response of your linearized model to the response of the nonlinear model using numerical integration to simulate the response of the nonlinear model.
  2. Compare the response of your linearized model to the real system's response by performing an experiment, and comparing data to the response of your linearized model.

In both options, you must take care to remember that your linearized model is a model for the dynamic response of the incremental output $\Delta y = y - \overline{y}$ as a result of applying an incremental input $\Delta u = u - \overline{u}$. Be sure to offset any linearized results by the NOP as appropriate!!

Model Validation: Numerical Integration

One way to deal with nonlinear system models, especially if all we want to do is visually see how they react to inputs, or compare a model with data, is to integrate their equations of motion numerically. Because the system is nonlinear, we can't use the linear differential equations techniques to easily obtain a closed-form solution to the system's behavior. Therefore, Numerical integration using something like Euler Integration or Runge Kutta Integration is usually a good option for investigating system response.

Given a general nonlinear vector-valued differential equation of the form $\frac{d\vec{x}}{dt} = f(\vec{x},t,\vec{u})$, where $\vec{u}$ are the system's input(s) and the system starts at some known initial condition, the Forward Euler solution involves solving for the function's value at discrete instants in time using the following equation:

\begin{equation} \vec{x}(t(k+1)) = \vec{x}(t(k)) + f(\vec{x},kT,\vec{u}(kT))T \end{equation}

You have probably seen Euler integration in more than one course so far. Perhaps you recognize it as "forward rectangular integration," which can be illustrated by the graphic below:

image-2.png

Note that in the equations describing Euler integration, computing $f(\vec{x},kT,\vec{u})$ is simply computing the value of derivative(s) at a particular time $t=kT$. Most times, a set of first-order differential equations describing a nonlinear system will manifest themselves as a list of scalar functions that depend on each other at that instant in time, along with the system's input(s) at that instant in time. Obtaining a solution to your system's governing differential equations using this method is easy, fast, and computationally 'cheap,' but does require an extremely small timestep $T$ to be accurate and numerically stable.

Using Numerical Integration to Simulate an Input-Output Differential Equation

To use Euler integration to solve a (possibly nonlinear) input-output differential equation of order $N$, you should first make sure that the input is defined so that there are no input derivatives, redefining your input as necessary. Then, you should solve your equation for the highest derivative. In other words, your equation should be written as:

$$\frac{d^N y}{dt^N} = f\left(\frac{d^N y}{dt^N},\frac{d^{N-1}y}{dt^{N-1}},\ldots,\frac{dy}{dt},y,u\right)$$

Then, Euler integration can be performed at every time $t=kT$ using a for-loop where $T$ is the integration time step, by using the differential equation to compute the highest derivative:

$$\frac{d^N y}{dt^N}(k-1) = f\left(\frac{d^N y}{dt^N}(k-1),\frac{d^{N-1}y}{dt^{N-1}}(k-1),\ldots,\frac{dy}{dt}(k-1),y(k-1),u(k-1)\right)$$

Then, also at every time $t=kT$, numerical integration can be performed successively on this highest derivative until the values of all derivatives of $y$, $\frac{d^{N-1}y}{dt^{N-1}(k)},\ldots,\frac{dy}{dt}(k),y(k)$, are estimated for each derivatives (including the 0th derivative, $y$) by computing the following for $n=N$ to $n=1$:

$$\frac{d^{n-1}y}{dt^{n-1}}(k)\approx \frac{d^{n-1}y}{dt^{n-1}}(k-1) + \frac{d^{n}y}{dt^{n}}(k-1)T$$

For example, $y(k)$, which could be written as $\frac{d^0y}{dt^0}$(k), would be estimated as:

$$y(k) = y(k-1) + \dot{y}(k-1)T$$

Example: Validating a Linearized Model against a Nonlinear Model

To illustrate how one might validate a linearized model against a nonlinear one using numerical integration, consider the linearization example we worked out earlier in the notebook. We will set up the simplest possible numerical integration code to evaluate a linearized model against a nonlinear one. Note that this is not your only option!! A dataset from an experiment can also be used to validate a linearized model!

In the code cell below, both the linearized model and the nonlinear model are simulated using Euler integration. The results of each simulation are compared in a plot for visual validation that the linearized model is a good approximation.

Because linearization requires a small deviation from the NOP, try modifying the code to use a larger $\Delta u$ value and watch how the agreement changes!

In [2]:
%Set up parameters of the models
K = 20;
a = 1;
%simulation parameters
simtime = 20;
dt = .001;
%simulation time vector-- we can make it common between both nonlinear and linear models
t = 0:dt:simtime;

%set up input situation
ubar = 10;
du = 0.5;
%set up the input vector for TOTAL input u = ubar + du
ubarvec = ubar*ones(size(t));
duvec = du*ones(size(t));
u = (ubar+du)*ones(size(t));

%set up variables for the nonlinear simulation
ynl = zeros(size(t)); %y for nonlinear
ydnl = zeros(size(t)); %ydot for nonlinear
%set up the initial conditions. Use the NOP equations to set them!
ybar = (K)/ubar;
ynl(1) = ybar;
ydnl(1) = 0;

%set up variables for the linear simulation. INCREMENTAL variables!!
Dyl = zeros(size(t)); % y for linear
Dydl = zeros(size(t)); %ydot for linear
%set up the initial conditions. Use the NOP equations to set them!
% For linearized, we are using Delta y, so the IC is 0!!
Dyl(1) = 0;
Dydl(1) = 0;


%perform Euler Integration for the NONLINEAR model
for k = 2:length(t)
    %highest derivative is a function of the others
    yddnl = 1/(a^2) * ( K - u(k-1)*ynl(k-1) - ydnl(k-1)* exp(-ydnl(k-1)) );
    %ydot can be integrated using ydd we just calculated
    ydnl(k) = ydnl(k-1) + dt*yddnl;
    %y can be integrated using the previous value of ydot.
    ynl(k) = ynl(k-1) + dt*ydnl(k-1);
end

%perform Euler Integration for the LINEARIZED model
for k = 2:length(t)
    %highest derivative is a function of the others
    Dyddl = 1/(a^2) * ( K - ubar*ybar - ubar*Dyl(k-1) - ybar*duvec(k-1) -Dydl(k-1) );%note Delta u (duvec) used here, not total u!
    %ydot can be integrated using ydd we just calculated
    Dydl(k) = Dydl(k-1) + dt*Dyddl;
    %y can be integrated using the previous value of ydot.
    Dyl(k) = Dyl(k-1) + dt*Dydl(k-1);
end

%shift linearized model up to start at NOP!
yl = Dyl + ybar;

% plot the results to see agreement!!
plot(t,ynl,t,yl)
xlabel('Time (s)')
ylabel('System Output (quarks)')
legend('Nonlinear Model','Linearized Model')

Note that as we attempt to simulate larger deviations from the NOP by setting $\Delta u$ to be a larger and larger step, the fidelity of the linearized model deteriorates!!

Exercise: Due by MIDNIGHT on Saturday, 9/25

Using the nonlinear version of the lumped-parameter, physics-based model at the beginning of this notebook, develop a linearized input-output differential equation model relating valve position. Assume that the value of the "nominal" valve resistance, or the resistance of the valve when valvepercent is just above 0% open, is $R_{nom} = 1.936x10^4 \frac{Pa \cdot s}{m^3}$.

Model Construction: Develop a Linearized Model Relating Valve Position to Turbine Speed

Show your model construction steps (linearization) in the markdown cell below. Your final product should be a linearized differential equation relating changes in valve position, $\Delta u_{v}$, to changes in turbine speed, $\Delta \Omega$.

You can use the following fixed assumptions:

  1. We will only model the system's behavior when the electrical portion of the system is disconnected.
  2. The damping between the turbine shaft and ground is aproximately 2,000 $\frac{Nms}{rad}$
  3. The turbine's rotational inertia $J$ is approximately 4250 $kgm^2$
  4. The turbine's displacement $D$ is approximately $1.0 \frac{m^3}{rad}$

YOUR ANSWER HERE

Model Validation

Using the simulator, collect a dataset that begins with the dam overflowing, and the valve open at 40% before data collection begins. When the user presses a button, the valve percentage should be increased to 50% while data are collected. Data should be collected for 5 seconds.

Controller Design for Data Collection Controller

Place a variable key and updated, properly-formatted FSM code ONLY below. Include comments to indicate what changes you made to your data collection code since the last assignment.

YOUR ANSWER HERE

Model Validation Using Data

Using your data for comparison, and either an analytical solution or numerical integration to simulate your linearized differential equation, validate your model by comparing it with your data on a single set of axes below.

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

error: No Answer Given!

Comment on the accuracy of your linearized model. How well does it work? If it works well, why do you think it does? If it does not work well, to what do you attribute the disagreement? Be analytical and use the tools we've developed so far in this course to answer this question.

YOUR ANSWER HERE