Table of Contents
- 1 ME480
- 2 Model Development steps
- 3 FSM
- 3.1 Controller Design: Boolean Algebra
- 3.2 Controller Design: The Four-Block Structure
- 3.3 Controller Design: The challenge of the first loop. How do you get in to the starting state?
- 3.4 Timers
- 3.5 Rising Edge Counters
- 4 Dynamic Modeling
- 5 Model Development: Constructing Models From Dynamic Tests
- 6 Model Development: Data-Driven System Models using Differential Equations
- 6.1 Model Order
- 6.2 Solving Differential Equations: the method of undetermined coefficients
- 6.3 Stability of Linear, Constant-Coefficient Differential Equations
- 6.4 Scoping a linear differential equation model using collected data
- 6.5 Constructing a model using collected data
- 6.6 Characterizing a response in terms of a differential equation model
- 6.7 First Order Linear Differential Equations
- 6.8 Characterizing First Order Systems
- 6.9 Second-Order Linear Differential Equations
- 6.9.1 Eigenvalues and characteristic equation
- 6.9.2 Free Response: general form for real, distinct characteristic roots
- 6.9.3 Free Response: general form for complex conjugate roots
- 6.9.4 Free Response: general form for repeated roots
- 6.9.5 Step Response with Zero Initial Conditions: Real, Distinct Roots
- 6.9.6 Step Response with Zero Initial Conditions: Complex Conjugate Roots
- 6.9.7 Step Response with Zero Initial Conditions: Repeated Roots
- 6.9.8 Step Responses with Nonzero Initial Conditions
- 6.10 Characterizing Second-Order Responses
- 7 Higher-Order Linear Differential Equations and the method of undetermined coefficients
- 8 Dynamic Physical Systems
- 9 Lumped Parameter Modeling of a Dynamic Physical System
- 10 Model Construction: Idealized, Lumped-Parameter Elements
- 10.1 The First Law of Thermodynamics
- 10.2 Power Form of the First Law of Thermodynamics
- 10.3 Power Variables: A-Type and T-Type
- 10.4 Idealized Sources
- 10.5 Idealized Lumped-Parameter Energy Storage Elements
- 10.6 Idealized Lumped-Parameter Energy Dissipation Elements
- 10.7 Idealized Lumped-Parameter Power Transducers
- 11 Model Construction: Interacting Idealized, Lumped Parameter Elements in a Physical System
- 11.1 Model Construction Tools: The Linear Graph Method
- 11.2 Elements in Linear Graphs: Sources
- 11.3 Elements in Linear Graphs: Dissipative and Energy Storage Elements
- 11.4 Elements in Linear Graphs: Idealized Transducers
- 11.5 Model Construction: Nodes
- 11.6 Model Construction: The principle of Continuity: Conservation at a Node
- 11.7 Model Construction: The Principle of Compatibility (Kirchoff's Voltage Law)
- 11.8 Model Construction: Building a Linear Graph Model of a System
- 12 Model Scoping Disciplined Process: Lumped-parameter Dynamic Physical Systems
- 13 Model Construction Disciplined Process: Lumped-parameter Dynamic Physical Systems
- 14 Model Construction: The Laplace Transform
- 15 Model Construction: The Transfer Function Representation of a System
- 16 Model Construction: Introduction to the Block Diagram
- 17 Disciplined Process: Block Diagram Simplification
- 18 Graphical interpretation of linear system behavior
- 19 Feedback Control
- 20 Anatomy of the closed-loop control system
- 21 The P(I)(D) Controller
- 22 Controller Design: P(I)(D) Controller Disciplined Process
- 23 Controller Design: Proportional (P) Control
- 24 Proportional-Integral (PI) Control
- 25 Common Controller Design Metrics for Closed Loop Systems Displaying Underdamped Second-Order Behavior
- 26 Proportional-Derivative (PD) Control
- 27 Controller Design: Root Locus Method
- 28 Summary of Root Locus Properties
- 29 Disciplined Process: Steps 1-5 of sketching the root locus
ME480¶
The ME480 Disciplined Process¶
You'll need to use a process to help guide decisions about what is needed and what might not be working to develop a controller for a system. The process can be extended from your system modeling experience in ME352. For both system modeling and controller design you need to develop an expectation of what will happen and then compare that to what does happen.
Model Development steps¶
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.
FSM¶
Finite State Machine (FSM) Model¶
However, to keep this model as simple as possible we could say "we assume the time it takes to transition between sleeping and awake is short enough that we aren't interested in the details of the transition". This means that we would only practically be interested in whether the brain was "sleeping" or "awake" as the transition would be so short it wouldn't be of interest.
This type of model is very common and allows you to represent different "situations" as Boolean States. A FINITE STATE MACHINE is a model that is always in a defined Boolean State and the time to transition between the states is neglected. While not everything can be accurately represented by a Finite State Machine model, they are very common in robotics and industrial control systems.
Each of the "Boolean States" in a finite state machine can take only one of two values:
- True (the machine is in this state)
- False (the machine is not in this state)
If you use a finite state machine to model your system then two assumptions are required
- The number of Boolean states is finite, meaning that you can count them.
- A Finite State Machine can be only be allowed to be in one Boolean state at a time.
Note that we are using the term "Boolean State" to distinguish from the "state" you're familiar with from system dynamics used in "state space" models.
FSM State Transition Diagrams¶
A state transition diagram is a design tool for FSMs that allows you to organize how your system will move from one state to another. It consists of bubbles representing each Boolean state, and arrows that show how the system is permitted to move between those Boolean states. A generic example is shown below-- note that it is not possible for the system to move from state 3 to state 1!
FSM - INPUTS, STIMULI (CONDITIONS) and TRANSITIONS¶
In general, it takes some kind of "stimulus" to cause a state machine to leave one Boolean state and go to the other. In our brain model the stimulus that brings you from "Sleeping" to "Awake" might be an alarm going off for your 8AM class. The transition from "Awake" to "Asleep" might be described by the condition that "you're tired and it's past your bedtime."
There are two types of stimuli that we will deal with in ME480:
- Inputs (stimuli from the user/outside world)
- Internal Stimuli (e.g. a timer internal to the controller)
Inputs are often things like switches, buttons, and control panels on a robot or machine. "Timers" in an FSM program keep track of how long it's been since some event has occurred, and "Counters" keep track of how many times something has happened. You'll learn about those in upcoming notebooks.
FSM: State Transition Tables¶
State transition tables are like a "key" that is used to read a state transition diagram. Each row represents a transition on the state transition diagram.
| Transition | Starting State | Condition | Ending State |
|---|---|---|---|
| A | Sleeping | Alarm rings | Awake |
| B | Awake | Bedtime and I'm tired | Sleeping |
Once this table is constructed, you can describe each entire transition from one state to another by reading across the table. For example, for transition B, we have:
"Transition B occurs when the state is 'Awake' and 'Bedtime and I'm tired' occurs."
Note that state transition tables are typically written using mathematical shorthand rather than full sentences. We will learn how to do this in a subsequent notebook.
FSM: OUTPUTS¶
In a FSM, what is considered an "output" is up to the engineer. Usually, an output is something that the machine uses to interact with its environment or its user... something observable to the outside world. Examples could be lights that indicate machine operation, signals that power motors or pumps, or text displayed on a user interface. One way to think about states is that they are "unique collections of system outputs."
What would the outputs be in our brain model?
The tools you need to design the gate-intersection system are same ones needed to control robots, AI characters in a video game, traffic lights, microwaves, digital watches, etc. Each of these machines and systems share the need to have a control system that makes a "decision" about how it will behave.
Synchronous FSMs¶
In this course, we will generally consider synchronous state machines, which basically means that their operation is governed by a clock signal executing one command after the other at a specified time. Generally, a synchronous state machine implemented in software consists of an "infinite" program loop that runs over and over. In each "loop" of the program, inputs are measured, logical conditions are evaluated to determine what state the machine should be in based on the user inputs and the machine's prior state. Then, at the end of each loop, outputs are written (think motors, lights, buzzers, etc.) based on what state the machine is in.
There are other types of state machines! For more information on types of state machine (this will not be on a test!) you can refer to this link or this one.
Because the state machines we will be implementing in this course are synchronous and running on an infinite loop, they need to make a decision about which state to be in every time through the loop. The decision has to be made regardless of whether anything of note has happened, and the software we write to implement a FSM must set each boolean state variable to either true or false on each pass through the loop.
Controller Design: Boolean Algebra¶
Boolean Variables: two possibilities¶
It's possible to intuit the concept of a Boolean variable. Here are some example manifestations:
| False | True |
|---|---|
| 0 | 1 |
| low | high |
| energized | de-energized |
| open | closed |
All of these words are equivalent. In some programming languages (like Arduino), you can use "0" or "false" or "LOW" interchangeably to represent "false," and "1," "true," or "HIGH" to represent true. Many languages, however, are more restrictive. MATLAB (or its open-source equivalent, Octave), for instance, insists that Boolean variables are assigned a value of either "true" or "false" or "1" or "0," respectively.
Operations on Boolean Variables¶
The building blocks of "Boolean Expressions" are the basic logical operators "OR" "AND" and "NOT." A Boolean expression can be thought of as a statement of equivalence comprised of Boolean (true/false) variables,
OR¶
The "or" operator is often written as a "+" sign. Example: "C is true if A or B is true" can be written:
\begin{equation} C = A+B \end{equation}In ME480 our simulators and programming languages use "||" to represent "or" so the expression above would be coded as:
C = A||B;
AND¶
The "and" operator is often written as a dot, $\cdot$, or ommitted entirely. Example: "C is true if A and B are true" can be written:
\begin{equation} C=A\cdot B \end{equation}In some circles (including this class), the $\cdot$ AND symbol can be omitted so the expression $Y=A\cdot B$ is the same as the expression $Y=AB$.
ME480 simulators and the programming languages used in the course use "&&" to represent "and" so the expression above would be coded as:
C = A&&B;
NOT¶
The "not" operator is often written as a bar over an existing variable, e.g. $\overline{A}$. For example, the statement: "C is true if A is not true" can be written:
\begin{equation} C=\overline{A} \end{equation}However, most programming languages do not allow you to use a $\overline{}$ sign to represent "not." In ME480, our simulators and hardware use languages that represent the "not" operator with a ! character. MATLAB (and Octave) are a bit different-- they represent "not" by either a tilde (~) or the word "not()." In an Arduino or ME480 simulator program, however, defining the Boolean variable "C" to be true if "A" is false, and false if A is true, we could write the following line of code:
ME480 simulators and Arduino use "!()" to represent "not" so the expression above would be coded as:
C = !(A);
MATLAB and Octive use "~" to represent "not" so the expression above would be coded as:
C = ~A;Common Symbols for Boolean Operations¶
The table below shows a summary of characters used for Boolean operators in different programming languages. It also includes the symbols use in other disciplines like philosophy.
Note that for "NOT," we include a "." to show that we are negating something. The "." is not part of the symbol(s).
| Language | OR | AND | NOT(.) | |
|---|---|---|---|---|
| ME480 | $+$ | $\cdot$ | $\overline{.}$ | |
| MATLAB | $ | $ | $\&$ | ~.,"not(.)" |
| Python | $|$, "and" | $\&\&$, "or" | ~., "not" . | |
| C,C++,Java,Arduino | $|$ | $\&\&$ | !. | |
| Philosophy, Logic | $\lor$ | $\land$ | $\neg .$ |
Order of Operations¶
Formally, the order of operations for Boolean Algebra is:
- NOT
- AND
- OR
Note: The "NOT" operator, denoted by the overline $\overline{.}$, has implied parentheses that span its length. Parentheses are used to group terms the same way as an overbar can be used, so parentheses and the "NOT" operator have equal precedence.
Much the same as it is in arithmetic and/or algebra. the "AND" ($\cdot$) operator takes precedent over the "OR" ($+$) operator, and parentheses can be used to group portions of an expression.
Truth Tables¶
When we evaluate a Boolean expression, it's often necessary to look at all possible combination of inputs (variables on the RHS of the expression) and evaluate their corresponding output (assigned variable on the LHS of the equation) to understand the implications of subjecting the expression to all possible combinations of inputs. To do this, we use a tool called the Truth Table. For example, the truth table for the expression $C=A\cdot B$ is given below:
| $$A$$ | $$B$$ | $$C=A\cdot B$$ |
|---|---|---|
| 0 | 1 | 0 |
| 1 | 0 | 0 |
| 0 | 0 | 0 |
| 1 | 1 | 1 |
Note that I needed to exhaust all possible combinations of the inputs in order to construct the table. As the number of inputs grows, so does the number of combinations you must test! In fact, there are $2^n$ combinations, where $n$ is the number of inputs in the expression!
Even with this complexity, truth tables are useful for testing more complex expressions "brute force." A good first example of this is the "XOR" or "exclusive OR" operator (shown as $\oplus$). This is like an "OR" operator but returns false if both A and B are true. By definition, the "XOR" operator is defined according to the following statement:
$$A\oplus B = \left(A+B \right)\cdot\overline{A\cdot B} $$The procedure for constructing a truth table is simple. Create a column for each variable in the Boolean expression, and then create columns for intermediate steps in computing the expression according to the order of operations. Make certain that you exhaust all possibilities for combinations of inputs, or else the truth table will not be "complete."
| $A$ | $B$ | $A \oplus B$ | $A+B$ | $AB$ | $\overline{AB}$ | $\left(A\right. + \left.B \right) \cdot \overline{A\cdot B}$ |
|---|---|---|---|---|---|---|
| 0 (false) | 0 (false) | 0 (false) | 0 (false) | 0 (false) | 1 (true) | 0 (false) |
| 1 (true) | 0 (false) | 1 (true) | 1 (true) | 0 (false) | 1 (true) | 1 (true) |
| 0 (false) | 1 (true) | 1 (true) | 1 (true) | 0 (false) | 1 (true) | 1 (true) |
| 1 (true) | 1 (true) | 0 (false) | 1 (true) | 1 (true) | 0 (false) | 0 (false) |
| ............... | ............... | ............... | ............... | ............... | ............... | ................................... |
The intermediate columns in the table above make evaluating the expression easier by grouping terms.
Properties of Boolean Expressions¶
It is possible to evaluate any Boolean expression using a truth table, but it is long and tedious! Using a few simplification rules presented below make the job becomes much easier.
Commutative Property¶
\begin{equation} A\cdot B = B\cdot A \end{equation}\begin{equation} A+B=B+A \end{equation}Associative Property¶
\begin{equation} \left(A\cdot B\right)\cdot C = A\cdot \left(B \cdot C\right) \end{equation}\begin{equation} \left(A+B\right)+C = A+\left(B+C\right) \end{equation}Distributive Property¶
The AND operator can be distributed into an OR expression (just like multiplication into addition!)
\begin{equation} A\cdot \left(B+C\right)=\left(A\cdot B\right)+\left(A\cdot C\right) \end{equation}Simplifying expressions: DeMorgan's Laws¶
DeMorgan's Laws, or the Laws of Negation, can be used to simplify long, complex Boolean expressions. It would be good practice to prove these using hand-written truth tables! Committing these to memory is actually quite usefule as they show up a lot in practical expressions common in logic-based programming and control (e.g. FSM design/implementation).
The first of DeMorgan's laws is one we like to call "neither!" It explains how to say that "neither A nor B can be true." \begin{equation} \overline{A+B}=\overline{A}\cdot \overline{B} \end{equation}
The second of DeMorgan's laws is one we like to call "Not this, or not that." It describes how to say that an expression is true if either A is false, or if B is false. Note that this is different than the case above, which says that neither A nor B can be true if the expression is to return true.
\begin{equation} \overline{A\cdot B} = \overline{A}+\overline{B} \end{equation}Example: Boolean Algebra Expression Simplification¶
Making use of these rules can help you simplify complex logic. Consider the following Boolean algebra expression.
\begin{equation} Y=A\left(\overline{A+B}\cdot B\right)+B \end{equation}The truth table for this would be epic. However, if we make the simplification for $\overline{A+B}$, we see that we get:
\begin{equation} Y=A\left(\overline{A}\cdot\overline{B}\cdot B\right)+B \end{equation}$B$ and $\overline{B}$ can never return true (B cannot be both true and false at the same time) so the whole parenthetical expression disappears, leaving us with $Y=B$.
Controller Design: The Four-Block Structure¶
One possible way to organize a FSM implementation (and the way that is required in ME480) is by separating the code into four blocks, which are given below. Using this structure makes it very straightforward to turn your state transition diagram into a functioning FSM implementation. The structure is given below.
Block 1: Process User Inputs & Sensors¶
Robots, factory assembly lines, CNC machines, handheld electronic games, HVAC systems, aircraft cockpits, and other automatic control systems whose operation can be supervised by a state machine usually have some way to interact with a human operator. Types of input devices are widely varied, but momentary switches and "latching" switches (toggle or locking switches like a light switch) are most common. In your laboratory assignment this week you have seen (or will see) that a switch can be used to generate a "high" voltage or "low" voltage. The outputs of these switches are the most basic type of Boolean input for a state machine design.
Block 1 is also where the program needs to take any "raw" inputs from the user and process/translate them for use in the state transition logic. For example, if your design requires you to hold a button for 2 seconds, you may need to write some logic that uses a "timer" (which we will cover later) to check how long a certain momentary button has been held.
Block 2: Evaluate State Transition Logic¶
This portion of the program is the most critical to implement correctly, but can be the easiest to implement if you have carefully considered the design of your state transition diagram.
To write the code for block 2, simply write a single line of code representing to the boolean expression for each transition in your state transition table without including the ending state. This means that in block 2, you should have as many lines of code as there are transitions in your state transition table. This is a simple but extremely powerful way to make a quick self-assessment of the reasonableness of your block 2 code.
Block 3: Set States¶
This section of code will finally "decide" which of our states will be true for this loop of the program. We look to the diagram to determine how to use the state transitions to set our current state. Specifically, we look at each state in our diagram and count the number of arrow heads pointing at it. Each arrow, as we know, represents one of the transition variables we defined in block 2. So all that is required now is to write the logic to say that a state is true if any of the state transitions ending in that state are true.
Functionally, this might manifest for a single state as:
State_x = T1||T2||T3;
where T1, T2, and T3 are the all state transition variables with an ending state of State_x. Let's write block 4 for our two-state program example.
Block 4: Activate Outputs For Current State and Update Variables¶
This (final) part of your Boolean Algebra program is fairly self-explanatory. Now that I know which state my machine is in, I'll use that information to activate any outputs (lights, sounds, motors, pumps) that are associated with this state. Because I'm "finished" with all of the program's tasks, I'll store the current state of relevant variables in an "OLD" variable so I will be able to access in the next loop of the program after inputs have been updated. In this case, we just need to store the "OLD" value of SW1 so that we can detect unique presses.
Controller Design: The challenge of the first loop. How do you get in to the starting state?¶
In a properly coded FSM, it is necessary to somehow initialize one of the design's states to true when the program first starts up. This is vital because all of the transitions in a state transition diagram depend on the context of the machine already being in a particular state.
Timers¶
In many applications where a FSM is a good choice for high-level program operation, you'll find yourself needing to design your state machines' transitions too use a few common stimuli. For instance: if you need to create a system that "waits" for a specified period of time after a button is pressed, or requires that a button is pressed for a particular length of time to "make sure" that the user meant to enter a particular machine state, you'll likely need to implement (or use) some kind of timer so that the stimulus of "wait for a specified amount of time" can be implemented in your state transition diagram and chart.
Timers are extremely common in finite state machine design. They behave in the following way:
- When the boolean input to the timer is TRUE, the timer counts up (usually in milliseconds).
- When the elapsed time is greater than some threshold, the timer's output variable is set to TRUE. Otherwise it is FALSE.
- When the input to the timer is FALSE, the timer stops counting and the elapsed time is reset to zero.
Rising Edge Counters¶
A rising edge counter helps a finite state machine keep track of how many times something has happened. This, like a timer, is a common need for designing the stimuli for transitions in FSMs. A standard counter that you might find on an industrial controller has three inputs: an "up input," a "down input," and a "reset input."
- If the user (or the state transition logic) sets the "up" input true with a unique rising edge (was false, now true-- like a button press), an internal variable (integer) in the counter will count up by one.
- If the using input sets the "down" input to be true with a "unique rising edge",it will bring the count (the counter's "accumulator") down by one. Generally, counters have a minimum value of 0.
- Whenever the "reset" input on the counter is true, the count is resets to zero.
- When the counter's accumulator reaches a preset value the output variable is set to TRUE. Otherwise it is FALSE.
Dynamic Modeling¶
Model Development: Mathematical Models¶
There are many strategies for developing equations of motion that describe the dynamics (motion) of a system. You've been exposed to several in your career as an engineering student. In this course, the dynamic models you will see are differential equations that relate a model's inputs to its outputs. In this course, we will use one of two general approaches to develop a system's governing differential equation:
- Data-driven approach: Conduct one or more experiments, and allow the experiments' results to guide model assumptions, complexity and structure. Validate the model using a different set of experiments that show the model has predictive power.
- Bottom-up approach: Begin with physics-based models of each component using appropriate assumptions. Build a model by connecting the system's components and capturing their interactions using the principles of continuity and compatibility.
Both of these approaches can be valid. Often, some combination of both strategies is necessary before a model is complete. Today, we will focus on the first steps towards building a model using a Data-driven approach.
Model Development: Constructing Models From Dynamic Tests¶
One of the most common ways to investigate a system's behavior when beginning to develop a model for its dynamics is to perform a test that involves changing the system's inputs and watching how its outputs evolve over time. By observing the relationships between the system's inputs and its outputs, we can often make determinations about what sort of mathematical model might be required to describe the system's behavior.
Step Response Tests¶
A "step input" is a common type of input used to investigate a system's dynamics. Colloquially, providing a system with a "step input" means that the system begins at an equillibrium defined by a steady-state input and a steady-state output. Then, the input is changed suddenly.
The unit step function (sometimes called the Heaviside step function after Oliver Heaviside) is defined mathematically as:
\begin{equation} u_s(t-t_0)=\left\{ \begin{matrix} 1 & \forall t\geq t_0 \\ 0&\forall t<t_0 \\ \end{matrix} \right. \end{equation}Applying a step change in a system's input is often represented mathematically by applying a scaled version of the heaviside unit step function, which by definition has a magnitude of 1. Generally speaking, a step function of magnitude $U$ that occurs at time $t_0$ can be written:
$$\Delta u(t) = U\cdot u_s(t-t_0)$$True Heaviside step functions represent instantaneous change in a system's input. Instantaneous changes are impossible in a physical system because they would require infinite power to achieve, but step functions are often a good approximation for "sudden" changes in a system's inputs. Step response tests are tests in which an approximate step change in a system's input is applied while the system's output(s) is/are measured.
Experimentally determining local system stability from a step response test¶
There are several mathematical definitions of "Stability" for a dynamic system. In this course, when we say that a system is stable, we generally mean that it is "Bounded Input Bounded Output (BIBO)" stable. What does this mean? It means that for any finite input, the system will produce a finite output.
Technically, "proving" that a system is globally stable in this way using experimental data would require giving the system every possible input and watching to make sure that its output stays finite. However, most systems are only subjected to a relatively small range of inputs. Consider your zumo or your lab rig's motor: you can only provide voltage inputs between 0 and 5V to the lab rig, so why worry about how the lab rig acts when it is fed 4,000 volts?
This is the concept of local stability, and without mathematical tools (which we will get to shortly), it's about all we can do to say that if a system is given step inputs with magnitudes within the range of interest, and the system reaches a steady state, it is likely stable. The figure below shows some examples of systems responding to step inputs:
Steady State in a step response test¶
For this purposes of this course, we define "steady state" behavior of a system in response to inputs that are nonperiodic (such as a step input) to occur when the output of the system is no longer changing. This definition changes slightly for periodic inputs (like sine waves).
Steady State Gain in a step response test¶
If a system reaches a steady state (if it is BIBO stable), a dynamic test can also give us information about the steady state ratio of output to input for the system. Think of this like a "calibration:" it tells us "how much output" we get "per unit input" for the system. Mathemtatically, we can define steady state gain for a step response test as the ratio of the change in output over the change in input for a step input as:
$$K_{ss} = \frac{y_{ss}-y_0}{u_{ss}-u_0} = \frac{\Delta y_{ss}}{\Delta u_{ss}}$$For other types of inputs, such as sinusoids, ramps, or impulse functions, the definition is slightly different, but the conceptual definition of steady state gain remains the same. It always tells us how much a system amplifies or attenuates its input.
Linearity: the principle of superposition¶
A system is linear if it satisfies the principle of superposition:
If an input change $\Delta u_1(t)$ produces an output change $\Delta y_1(t)$ and an input change $\Delta u_2(t)$ produces an output change $\Delta y_2(t)$, then
input change $\Delta u(t) = c_1\Delta u_1(t) + c_2 \Delta u_2(t)$ produces an output change $\Delta y(t) = c_1\Delta y_1(t)+c_2\Delta y_2(t)$ for all pairs of input changes $\Delta u_1(t),\Delta u_2(t)$ and constants $c_1,c_2$.
This means that the system satisfies the principle of superposition.
One key consequence of the principle stated above is that if one doubles the input change $\Delta u(t)$ for a system, a linear system will produce exactly double the output change $\Delta y(t)$ at every time $t$.
Unfortunately, essentially all physical systems known to humans are both nonlinear and of essentially infinite order. However, given that we're often only concerned with the behavior of a system within certain small regions of inputs and initial conditions, many physical systems are approximately linear in a particular region of interest to an engineer.
2% Settling Time¶
Settling time gives us an idea of how long the transients (changing behavior) of a dynamic system last. We use settling time to classify a system's dynamics as "fast" or "slow" relative to our control system's design goals, or when comparing one configuration of a system to another.
We will link settling time to mathematical models and their properties, but it is an empirical concept that can be obtained simply by analyzing a dataset. The 2% settling time is defined as the last time at which the system's output change $\Delta y(t) = y(t) - y_0$ has an absolute value that is greater than 2% away from $\Delta y_{ss} = y_{ss}-y_0$. This is illustrated graphically below for two systems, one of which is oscillatory, and another which reaches its steady state value asymptotically.
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.
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:
- 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.
- Find the particular solution $y_p(t)$ to the differential equation, where the system is subjected to a particular input ($u(t)\neq0$)
- Find the total response of the system by summing the homogeneous and particular solutions.
- 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:
- 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?
- Collecting and plotting data of the output (and input, if possible) during the test
- Confirming that the system in question appears to be stable and is approximately linear for ranges of inputs you're interested in
- 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:
- Finding the approximate eigenvalues your differential equation model will need in order to match your system's transient response.
- 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¶
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¶
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)$.
Settling Time¶
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)$$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¶
Free Response: general form for complex conjugate roots¶
Free Response: general form for repeated roots¶
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¶
Step Response with Zero Initial Conditions: Complex Conjugate Roots¶
Step Response with Zero Initial Conditions: Repeated Roots¶
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.
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.
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:
- Underdamped, with a complex conjugate pair for the eigenvalues $p_1,p_2$
- Critically damped, with two repeated real eigenvalues $p_1 = p_2$
- 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$.
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:
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}$$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.
Dynamic Physical Systems¶
A dynamic physical system is a collection of interconnected, physical components. At least one of these components must store energy for the system to display dynamic (time-varying) behavior in the absence of a time-varying input.
Lumped Parameter Modeling of a Dynamic Physical System¶
Lumped parameter modeling is the act of approximating a real, physical system comprised of interconnected real, physical components, each of which may have spatially or temporally varied physical properties as a system of interconnected idealized components that each have one "lumped," spatially-and-time-invariant physical property.
This is the approach to physics-based modeling we will use in ME480.
Model Construction: Idealized, Lumped-Parameter Elements¶
The First Law of Thermodynamics¶
When we go to construct dynamic, physics-based models for lumped-parameter, idealized elements (which are themselves one-element systems), we can start building our model by thinking about how the element stores and/or dissipates energy when work is applied to it.
You have seen the law of conservation of energy for a system, also called "the first law of thermodynamics," in your thermodynamics course. It states that a change in a system's internal or stored energy from energy state "1" to energy state "2" $E_{1\rightarrow 2}$ must be caused by either heat transfer $Q_{1\rightarrow 2}$ into our out of the system boundaries, or by work $W_{1\rightarrow 2}$ done to or by the system.
$$E_{1\rightarrow 2} = Q_{1\rightarrow 2} - W_{1\rightarrow 2}$$For this discussion, we will use the "Mechanical Engineering" convention for signs, in which heat transfer into the system is defined positive and work done by the system is defined positive. Note that this means that if a system is doing "positive work" on its surroundings, its stored energy will decrease in the absence of heat transfer. If positive heat transfer occurs at the system boundaries, the system's stored energy will increase. Other sign conventions are possible, and they vary from field to field.
The form of stored energy $E$ and the work done to or by the system $W$ will vary based on whether we are discussing a mechanical, electrical, fluid, thermal, or mixed system or element, but as you are probably aware, the units of all of these types of energy are dimensionally consistent with $\frac{kgm^2}{s^2}$ or "Joules."
When we approximate the physical construction of a system using idealized, lumped-parameter elements, we often say that these elements fall into one of the following categories:
- Idealized Sources: Deliver power to the system by providing a known value of one of the two variables in the system's power equation, regardless of how much power is required.
- Energy Storage Elements: store energy without heat transfer
- Dissipative Elements: dissipate energy without storing any
- Power Tranducers: transform the energy's type without storage or heat transfer
All of the commonly used lumped parameter elements we will discuss fall into one of these three categories. If a system or component does more than one of the three things mentioned above, it can often be split up into idealized elements that only perform one job.
Power Form of the First Law of Thermodynamics¶
Because we are working towards building physics-based dynamic models of a system, and are interested in how our system's behavior evolves over time rather than just between two energy states "1" and "2," we can shrink the duration over which our first law equation is applied, and look at the rate of energy change in smaller and smaller time intervals. In the limit of $\Delta t \rightarrow 0$, we end up taking the derivative of the first law of thermodynamics. Because the derivative of energy is power, with SI units of Joules/second, we call this the power form of the first law, and it can be applied to a system at any arbitrary moment in time:
$$\dot{E} = \dot{Q} - \dot{W}$$Power comes in many forms-- in this course, we will focus mainly on electrical, mechanical, and incompressible fluid power (flow work), which are the products of two key variables in each case:
- Electrical Power $\dot{E}_{elec} = V\cdot i$, a product of Voltage $V$ and current $i$
- Translational Mechanical Power $\dot{E}_{mech,t} = v\cdot F$ with velocity $v$ and force $F$
- Rotational Mechanical Power $\dot{E}_{mech,r} = \Omega \cdot T$, with angular velocity $\Omega$ and torque $T$
- Incompressible Fluid Power $\dot{E}_{fluid} = P \cdot \dot{\mathcal{V}}$, with pressure $P$ and volumetric flow rate $\dot{\mathcal{V}}$.
Power Variables: A-Type and T-Type¶
In each of the cases above, the variables in the power equation can be grouped into categories based on how they are measured. these two categories are "across" and "through."
- Across-type (A-type) variables must be measured relative to a reference. Only differences in these quantities between two points can be measured. These variables are pressure $P$, velocity $v$, angular velocity $\Omega$, and voltage $V$.
- Through-type (T-type) variables must be measured by placing an instrument in series with a system element. Examples include current sensors, load cells, and flow rate gauges. The T-type variables in incompressible fluid, mechanical, and electrical systems are volumetric flow rate $\dot{\mathcal{V}}$, force $F$, torque $T$, and current $i$.
Idealized Sources¶
It is sometimes safe to assume that the input to a dynamic, physical system is "ideal" in an energetic sense. Idealized sources provide power to the system, and can change its energetic sense.
Assumptions¶
Idealized sources can provide a known input to a system, which is usually a power variable, regardless of how much power is required to maintain that known input. For example, if the input to a system is force, then we might say that an "idealized force source" is able to provide a as much power as is required to the system to maintain a known input force. These infinite power sources are not real, but many times are a good approximation for systems that operate in a relatively small range of energetic states.
Because power is the product of a T-type and an A-type variable, idealized sources are often classified as either "T-type" sources, which provide a known T-type input regardless of what A-type variable is required, or "A-type" sources, which provide a known A-type variable regardless of what T-type variable is required.
Note that in the figure, a line representing the possible behavior of a "real" power source, which only has a finite available power, is also included. Many real power sources can be approximated as ideal T-type or A-type sources if the power required by the system is low relative to the power source's capability.
List of Common Idealized Power Sources¶
Common idealized T-type sources include:
- current sources in electrical systems
- force or torque sources in mechanical systems
- volumetric flow rate sources in fluid systems
Common idealized A-type sources include:
- voltage sources in electrical systems
- velocity or angular velocity sources in mechanical systems
- Pressure sources in fluid systems
Idealized Lumped-Parameter Energy Storage Elements¶
$$\require{cancel}$$Idealized energy-storage elements are a useful approximation of many real objects that store energy. Stored energy comes in many forms: mechanical kinetic and potential energy, chemical energy, and thermal energy are all examples. If an object primarily stores only one type of energy, it may be a good candidate to treat as an idealized energy storage element.
Assumptions¶
Idealized energy storage elements (except those that store thermal energy explicitly) are assumed to have no heat transfer in or out. The net work done on or by the element must balance with its stored energy, which can be written formally using the first law of thermodynamics by ignoring heat transfer.
$$\dot{E} = \cancel{\dot{Q}} - \dot{W}$$Further, the lumped-parameter, idealized energy storage elements used to construct dynamic physical models are assumed to store only one type of energy, and to exchange energy with their surroundings using only one type of work.
For each particular type of energy storage element, these assumptions have different consequences. For fluid capacitors, any kinetic energy due to fluid entering or leaving the capacitor is ignored. For springs, the mass of the spring is ignored. For idealized rotational inertias, any material elasticity that could store potential energy is ignored. The list goes on, but the key thing to remember is that while no real object is actually an idealized, lumped-parameter energy storage element, many real objects are well-approximated by these assumptions because one form of energy storage vastly dominates any other relevant terms in the first law equation above.
Energy Storage Element Types¶
Energy storage in fluid, mechanical, and electrical systems is usually accomplished by accumulating either the T-type or the A-type variable in the power equation (not both). In mechanical systems, a spring stores potential energy as $E = \frac{1}{2}K x_{12}^2$, where $x_{12}$ is the spring deflection. Substituting Hooke's law into this equation yields $E = \frac{1}{2K} F ^2$, where $F$ is the force in the spring. Because the energy equation can be written in terms of the T-type variable $F$, the spring is considered a T-Type energy storage element. Conversely, an electrical capacitor, which stores energy as $E = \frac{1}{2} C V_{12}^2$, is considered an A-type energy storage element because voltage is an electrical system's across-type power variable.
Classifying energy storage elements this way allows us to draw analogies between different system types. It allows us to treat mechanical springs similarly to fluid inertors, capacitors similarly to masses, and so on.
Elemental Equations for Idealized Energy Storage Elements¶
An idealized energy storage element's "elemental equation" is a restatement of the first law of thermodynamics in power form. An elemental equation uses an empirical relationship, e.g. Hooke's Law in the case of a mechanical spring, and combines that relationship with the first law to tell us how the element's energy is accumulated in terms of the element's lumped parameter, or dominant physical characteristic.
List of Idealized Lumped-Parameter Elemental Equations¶
A list of the elemental equations for the lumped-parameter energy storage elements we may encounter in ME480 is shown below, along with the equation for the element's stored energy.
Idealized Lumped-Parameter Energy Dissipation Elements¶
$$\require{cancel}$$Idealized energy-dissipation or "dissipative" elements are a useful approximation of many real objects that do not store significant energy, but for which the net work done on or by the object is not zero. In order to satisfy the conservation of energy, this type of idealized element must transfer energy as heat to the environment, which is why they are called dissipative elements. They result in energy leaving the system boundaries as heat.
All real systems dissipate energy. If they did not, we would have perpetual motion machines! Idealized, lumped-parameter dissipative elements are often used to model the major dissipative processes and components in real systems.
Assumptions¶
Idealized dissipative elements are assumed to store no energy. The net work done on or by the element must balance with the heat transfer at the element's boundary. This can be written formally using the first law of thermodynamics by ignoring all energy storage.
$$\cancel{\dot{E}} = \dot{Q} - \dot{W}$$Further, the lumped-parameter, idealized energy dissipation elements used to construct dynamic physical models are assumed to exchange only one type of work with their surroundings. This could be electrical work, flow work, etc., but not a combination of these types.
Elemental Equations for Idealized Energy Dissipation Elements¶
An idealized energy dissipation element's "elemental equation" is a restatement of the first law of thermodynamics in power form. An elemental equation uses an empirical relationship, e.g. Ohm's law in the case of an electrical resistor, and combines that relationship with the first law to tell us how the element transfers heat outside of the system boundary in terms of the element's lumped parameter, or dominant physical characteristic.
List of Idealized Lumped-Parameter Elemental Equations¶
A list of the elemental equations for the lumped-parameter energy dissipation elements we may encounter in ME480 is shown below, along with the net power consumed by each element, which all must be transferred to the element's surroundings as heat.
Idealized Lumped-Parameter Power Transducers¶
$$\require{cancel}$$Idealized, lumped-parameter power converting transducers are often used to represent physical objects in a system that transform energy from one form to another. Motors convert electrical work to mechanical work. Pumps convert fluid flow work into mechanical work. Gears convert mechanical work at one angular velocity to mechanical work at another angular velocity. The key characteristics of power-converting transducers are that the input work and the output work are the same, meaning that the idealized transducer neither stores nor dissipates energy.
Assumptions¶
The power-converting transducer only does "one job," and that is power conversion. Therefore, if something can be said to be adequately modeled as a power-converting transducer, it cannot store energy or transfer it to the system's surroundings. In other words, its first law of thermodynamics equation in power form looks like this:
$$\cancel{\dot{E}} = \cancel{\dot{Q}} - \dot{W}$$These simplifications mean that the net work on the transducer must be zero-- in other words, the power into the transducer and the power out of the transducer must be the same.
$$\dot{W}_{in} = \dot{W}_{out}$$For a gear train, the assumption of no energy storage would mean that the gears must be massless. The assumption of no heat transfer would mean that the gears have no friction or damping in their bearings. For a motor, the assumption of no energy storage would mean that the motor shaft has no inertia. The assumption of no heat transfer would mean that the motor's armature has no electrical resistance, and that its rotating assembly has no damping. Are these reasonable assumptions?
Probably not. Approximating real transducers using only an idealized transducer element is often a mistake-- most real motors, pumps, and gears have losses and/or intrinsic inertias that store energy. These "extra" pieces of a real power transducer can be represented by breaking the real power transducer into a couple of lumped-parameter idealized elements. Often, the full set of elements needed to fully describe a real power transducer include energy storage and/or dissipative elements along with an idealized transducer. This separates each of the real energy storage, conversion, and dissipative processes in the transducer into "chunks."
List of Idealized Lumped-Parameter Elemental Equations¶
Model Construction: Interacting Idealized, Lumped Parameter Elements in a Physical System¶
Model Construction Tools: The Linear Graph Method¶
You will almost never scope a physics-based model for a system using only one idealized lumped parameter element. However, as the number of elements in your model increases, it can be increasingly difficult to keep track of how the elements influence one another and transfer energy to each other or to the environment.
One convenient and powerful way to help visualize how power (and thus energy) flow in a network (system) of interconnected lumped-parameter elements is to represent a system as a kind of "circuit," regardless of whether the system in question is electrical, mechanical, fluid, or mixed.
This creates a problem, in that it is hard to draw masses, tanks, or viscous dampers as "circuits" in the traditional sense. To make this operationally easier, system dynamicists have come up with several systems for drawing generalized circuit-like networks for use in understanding how a system's components interact. We will explore such a "universal" method to help operationalize model construction in ME480.
The approach we will use is based on prof. Seeler's system dynamics textbook. The methodology is called the Linear Graph method, and it was developed at MIT in order to treat all systems of one-dimensional, idealized, lumped-parameter elements in a consistent way. The law of the land in the linear graph method is "no special symbols," so while everything is represented as a circuit-like network of elements connected by "nodes" (points of constant across-type variable), the symbol for an electrical resistor looks the same as a symbol for a mechanical damper.
Elements in Linear Graphs: Sources¶
In a linear graph, idealized sources are represented as lines coming from "ground," which represents the zero reference point for the Across-type variable relevant to the type of system in question (velocity, pressure, or voltage). The source can be either an idealized Through (T) type source or an idealized Across (A) type source.
Elements in Linear Graphs: Dissipative and Energy Storage Elements¶
Lumped parameter energy storage and energy dissipating elements in linear graphs are drawn as simple arrows. These elements can all be thought of as having an "input port" and an "output port" for the through-type variable. Why?
Think about a spring or a damper, which both physically have two connections, one at each end. A section analysis and free body diagram of an idealized, massless spring or damper would show that the net force is the same throughout the element, so force is thought of as flowing through the element from one connection, or "node," to the next.
Any two-port element in a linear graph is drawn this way. Some examples are shown below.
The only time a lumped-parameter "two port" element is drawn differently is in the case of an A-type energy storage element that really only has one physical connection. For example, fluid capacitors store energy in pressure that is always measured with respect to the surrounding pressure, which is the "ground node" or zero reference point for the A-type variable in the system. Similarly, translational masses and rotational inertias store energy in velocity, but this velocity must be referenced to an inertial frame in order for Newton's laws to apply. Therefore, they too must always be connected to "ground." The tricky part is that these elements aren't physically connected to the ground reference, so we draw them with a dotted line to show that they are referenced to ground but not physically connected. Examples are shown below.
Elements in Linear Graphs: Idealized Transducers¶
In a linear graph, idealized transducers connect portions of the diagram that represent different power types-- for example if we are talking about a motor, the power into the motor has an across type variable of voltage and a through type variable of current. The power out of the motor has an across type variable of angular velocity and a through type variable of torque. The generic symbol for an idealized transducer is shown below:
Note that the "ground" reference is not connected between the two "sides" of the transducer. This is to symbolize that the A and T type variables on each "side" of the transducer might have different units. It's also worth mentioning that a transducer is the only type of element we will use in a linear graph that has four "ports" rather than two-- it has two ports for its "input" power, and two ports for its "output" power.
Model Construction: Nodes¶
In system dynamics, a "node" is a fictitious element in your system through which power is transferred without energy storage or "losses" (energy transfer to outside the system boundary). In electrical, fluid, and mechanical systems, a "node" represents a fictitious place where no heat transfer occurs from the system to its environment and no energy is stored.
In terms of the first law of thermodynamics, these assumptions result in the following:
$$\require{cancel}$$$$\cancel{\dot{E}} = \cancel{\dot{Q}} - \dot{W}$$This means that the net work on the node must be zero.
We often imagine that idealized, lumped-parameter elements in our system model are connected to one another through this type of lossless interface. the node "splits" the through-type variable in the power equation, distributing it to the elements to which it is connected. A single node will by definition be a place in your system where the across type variable in the power equation is constant. Jumper wires and breadboard rows are examples of "nodes" in electrical circuits. Rigid connections between a mass and a spring could be conceived of as "nodes" in a mechanical system. A short section of pipe with negligible fluid resistance where multiple flows split or come together in a fluid system also might be well approximated by a "node" of constant pressure.
Model Construction: The principle of Continuity: Conservation at a Node¶
The principle of continuity is a re-statement of the first law of thermodynamics with specific assumptions. For many types of systems, the assumptions boil the first law of thermodynamics down to a simple statement about conservation of mass. In these cases, the principle of continuity states:
"flows into a node must balance with flows out of a node,"
Where the definition of "flows" depends on the type of system at hand. For electrical systems, the principle of continuity at a node boils down to a "conservation of current," or Kirchoff's Current Law:
$$\sum i _{node} = 0$$With currents into a node considered positive, and currents out of a node considered negative.
For a mechanical system, Force is the quantity that "flows through" elements via Newton's 3rd law (making an imaginary cut anywhere along a rigid body and drawing a free body diagram will confirm this). Thus, a "node" in a mechanical system is a rigid, imaginary element with "no" mass. Often, a node is superimposed on top of a translating or rotating idealized mass, and we imagine that the inertial reaction force from that idealized mass flows into or out of the node as appropriate. This is similar to the use of the D'Alembert principle which conceives an "inertial reaction force" exerted by a mass rather than thinking about the familiar form of Newton's second law, $\sum F = m\dot{v}$. The D'Alembert vs. Newtonian viewpoint is probably responsible for the war against "centrifugal force" waged in many high school physics classrooms across the U.S.
Writing the first law of thermodynamics for a node in a mechanical system, and canceling that its (single) velocty from the first law equation, one obtains:
$$\sum {F}_{node} = 0$$For a fluid system (incompressible), a "node" might represent an infinitesimal junction between two pipes or valves, or a junction between a valve and a tank. Writing the first law of thermodynamics for this "infinitely thin" element in which only one pressure exists (and can thus be canceled), one obtains:
$$\sum \dot{\mathcal{V}}_{node} = 0$$With volumetric flows into the node being considered positive, and flows out of the node negative.
Model Construction: The Principle of Compatibility (Kirchoff's Voltage Law)¶
The principle of continuity states that the sum of the voltage drops around any closed "loop" in a circuit-like network must be zero. This is a direct consequence of the first law of thermodynamics. Consider the following circuit, for which 3 closed loops exist. The circuit diagram and the linear graph are shown side by side for comparison. They contain exactly the same information.
For Loop 3, the principle of compatibility states: $$\require{cancel}$$ $$ V_{g1} + V_{12} + V_{23} + V_{34} = 0$$
Which, given the definition of $V_a - V_b = V_{ab}$, could be re-written as:
$$ \cancel{V}_{g} - V_1 + V_1 - V_2 + V_2 - V_3 + V_3 - \cancel{V}_{g} = 0$$Similarly, for Loop 1, the principle of compatibility states:
$$V_{g1} + V_{12} + V_{2g} = 0$$And for Loop 2, the principle of compatibility states:
$$V_{g2} + V_{23} + V_{3g} = 0$$Compatibility can also be appled to fluid systems (replacing voltages with pressures) and to mechanical systems (replacing voltages with velocities).
Model Construction: Building a Linear Graph Model of a System¶
The aim of the linear graph method is to take a diagram showing physical components of a system and to turn it into a circuit-like network of idealized, lumped-parameter elements. The model must be properly scoped before a linear graph model is constructed. This means that the model's inputs, outputs, and its list of constitutive idealized elements must be known before a linear graph is constructed.
To construct a linear graph, you can follow the general procedure below.
- Draw a ground reference for each type of power in the system (rotational mechanical, translational mechanical, electrical, or fluid)
- Draw and label all of the nodes in the system. Look for unique nodes by looking for places in the system with a unique value of the A-type power variable for each respective element (velocities, angular velocities, pressures, and voltages).
- Connect the nodes using your idealized, lumped parameter elements, including transducers. Ensure that each end of each element has the correct A-type variable value.
- Annotate the model with unique values of the T-type power variables in the system (currents, forces, torques, flow rates). Assume directions for these T-type variables according to your best guess of the direction of power flow from input to output.
Once your linear graph representation of your system is complete, you are ready to continue with model construction. The linear graph can be treated like an electrical circuit. The principles of continuity (KCL) and compatibility (KVL) will be helpful in building a differential equation or set of differential equations describing your system.
Model Scoping Disciplined Process: Lumped-parameter Dynamic Physical Systems¶
Model scoping for a dynamic, lumped-parameter model can start similarly to model scoping for a purely empirical model.
- Begin by choosing an input and output for the system. Determine whether an idealized A or T type input is most appropriate. The units of the system's input and output will often tell you if power transduction must occur inside the system. For instance, if an input is electricity, and an output is fluid flow, these represent two different types of power, so conversion must take place inside the system.
- If possible, a step response, initial condition response, or other dynamic test should be performed on the system to get an idea of the system's order.
- The apparent order of the system will be equal to the number of independent energy-storing elements that are required for inclusion in the system. Two energy-storage elements are independent if their energy storage equations cannot be written in terms of the same across or through-type variable.
- Look at the physical components in the system. Determine which of these components are likely to store significan energy, and represent the energy storing processes in your system as idealized energy storage elements. In mechanical systems, look for masses and inertias. In electrical systems, look for components that may have significant capacitance or inductance. In fluid systems, look for reservoirs or tanks to model as capacitors, and consider modeling long pipes with significant fluid mass inside them as fluid inertors. Start simply and allow yourself to be guided by your conclusions about model order above.
- Look at the physical components in the system. Determine where there are likely to be dissipative processes in the system, and represent these as idealized dissipative elements. In mechanical systems, look for damping or friction between two elements, or between an element and ground. In electrical systems, look for resistors or elements that may have resistance. In fluid systems, look for valves, long pipes, or orifices that may significantly resist flow.
- Guided by step 1, determine whether any power transduction takes place. Choose idealized power transducers to include in your system as appropriate.
Model Construction Disciplined Process: Lumped-parameter Dynamic Physical Systems¶
After you have scoped your lumped-parameter dynamic model, it is time to move on to model construction. What follows is one possible process you could follow. There are others, but this is the one we will refer to in ME480. You have followed a similar process in ME352, so we will not provide exhaustive, detailed examples here.
- Describe the behavior of each element in your model (write its elemental equation and energy storage equation if appropriate)
- Determine if the number of independent energy storing elements in your model matches your expectations from model scoping. If it does not, you may want to reconsider your scoping choices.
- Draw your system of interconnected, lumped-parameter, idealized elements as a circuit-like network using a tool like the Linear Graph method. This will show how each idealized, lumped-parameter element in your model is connected by "nodes" of constant across-type variable.
- Write the equations of continuity for each node in your model.
- Write the equations of compatibility for each independent "loop" in your model.
- Using algebra, treat the set of elemental, node, and loop equations as a system of coupled equations. Reduce this set of equations to a suitable form so that they can be solved and/or simulated. "A suitable form" is your choice, but may be an input-output differential equation, a transfer function, or a set of coupled differential equations in "state space" form.
- Confirm that your final model (in the form of a differential equation) has the correct order, is stable for your elements' parameter values, and that the units of each term in your differential equation(s) are consistent.
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:
- Solve (or simulate) the nonlinear differential equation.
- Replace any nonlinear physical elements with approximately equivalent linear ones.
- 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:
- 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.
- Using the Taylor Series expansion, linearize all non-constant, nonlinear terms in your system's differential equation.
- Substitute the linearized terms back into your differential equation, substituting all remaining $y = \overline{y} + \Delta y$ and $u = \overline{u} + \Delta u$.
- 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.
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:
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$$Model Construction: The Laplace Transform¶
The Laplace transform is an analytical tool that allows us to turn linear, differential equations into algrebraic equations. This makes them easier to modify, combine, and manipulate. The Laplace transform also allows us to use convenient tables to obtain solutions for differential equations. The laplace transform of a function in the time domain is defined as:
$$F(s) = \mathcal{L}\left|f(t)\right| = \int_0^{\infty} f(t) e^{-st}dt$$Taking the Laplace transform of a differential equation turns it from an equation with time derivatives into an equation that is a polynomial in "s." For example, assuming zero initial conditions, the term $\frac{d^n y(t)}{dt^n}$ would turn into the term $s^n Y(s)$. In general, including the effects of initial conditions, we can write:
$$\mathcal{L}\left| \frac{d^n y(t)}{dt^n} \right| = s^n Y(s) -s^{n-1} y(0) - s^{n-2}\frac{dy}{dt}(0) - \ldots - s\frac{d^{n-2}y}{dt^{n-2}}(0) - \frac{d^{n-1}y}{dt^{n-1}}(0)$$Model Construction: The Transfer Function Representation of a System¶
In this course, we'll often need to combine systems to create new ones, or manipulate system models during a controller design process. For this reason, we usually will represent our system as a transfer function, assuming that the system is linear, it has only one scalar input and one scalar output. A transfer function is simply the ratio of a system's output/input.
To obtain the transfer function for a system, assume that its initial conditions are zero. Then, take the Laplace transform of your system's input-output differential equation, and solve for the ratio $\frac{Y(s)}{U(s)}$ where $Y(s)$ represents the system's output, and $U(s)$ represents the system's input.
This procedure will result in a rational function of two polynomials. For a differential equation with $M$ input derivatives and $N$ output derivatives, this looks like:
$$\frac{Y(s)}{U(s)} = \frac{b_m s^M + b_{M-1} s^{M-1} + \cdots + s b_1 + b_0}{a_m s^N + a_{N-1} s^{N-1} + \cdots + s a_1 + a_0}$$This rational expression contains all of the information we need about our dynamic system model. It allows us to look at the model as an "operator" that can produce output as its function of the input given to the system.
Zeros of the system are the roots of the numerator polynomial, and dictate how the system responds to changes in input.
Poles (eigenvalues) of the system are the roots of the denominator's polynomial, and tell us everything we need to know about the system's intrinsic transient characteristics.
Tools for analyzing transfer function behavior¶
Finding the zeros of a transfer function¶
The zeros of a transfer function can be found by setting the numerator of the transfer function equal to zero. Zeros do not affect stability. They merely show us how a system reacts to sudden changes in inputs. Consider the following transfer function:
$$\frac{Y(s)}{U(s)} = \frac{s+b_0}{s+a_0}$$by the rule above, we can see that this transfer function has one zero at $s = -b_0$. but what does this zero tell us? To understand it, let's simply break the transfer function into two pieces:
$$\frac{Y(s)}{U(s)} = \frac{s}{s+a_0} + \frac{b_0}{s+a_0}$$So our system can be conceptualized as consisting of two separate systems. By cross-multiplying to study the expression we get:$$Y(s) = \frac{1}{s+a_0} \left(s U(s) + b_0 U(s)\right)$$ Here one "piece" of the system reacts to the input itself $U(s)$, and the other reacts to the derivative of the input, $s U(s)$.
By replacing $s$ with $-b_0$, the larger the magnitude of the zero $\left|-b_0\right|$, the less of a contribution the derivative piece, $sU(s)$ will have on my system's overall response. The smaller the zero's magnitude is, the more important the derivative of the input becomes!!
Initial Value Theorem: Find a model's initial response to an input¶
Another nice thing about a Laplace Transform representation of a system, especially if that system is in transfer function form, is that there are quick, easy ways to get information about our system's response to an input. The first thing we'd like to know is: at $t=0^+$, what is the value of the system? It started at rest, but it's possible that our input "instantaneously" changed something about the system. This is what the initial value theorem is great for. Let's assume we have a system transfer function $\frac{Y(s)}{U(s)} = T(s)$. Then, the response of our system to a particular input is $Y(s) = U(s)T(s)$.
To find the system's response to the input at time $t=0+$, or just after the input has been applied, we can use the Initial Value Theorem as long as our system is BIBO stable.
If the Laplace transform of a $y(t)$ is $Y(s)$, then the initial value theorem states:
\begin{equation} \lim_{t\rightarrow 0} y(t) = \lim_{s\rightarrow \infty} sY(s) \end{equation}The Final Value Theorem: Find a model's steady-state, final value¶
Just like the initial value theorem, the final value theorem allows us to obtain the system's steady-state value given its laplace transform representation. Let's assume we have a system transfer function $\frac{Y(s)}{U(s)} = T(s)$. Then, the response of our system to a particular input is $Y(s) = U(s)T(s)$.
To find the system's response to the input at time $t=\infty$, we can use the Final Value Theorem as long as our system is BIBO stable.
If the Laplace transform of $y(t)$ is $Y(s)$, and if $sY(s)$ is analytic on the imaginary axis and in the right half of the $s$-plane, then
\begin{equation} \lim_{t\rightarrow \infty} y(t) = y_{ss} = \lim_{s\rightarrow 0} sY(s) \end{equation}The final value theorem is particularly useful for finding a system's steady state gain in response to a particular input (such as a step, although this is not the only input for which a steady state gain may exist). The steady state gain, as before is defined as $K_{ss} = \frac{\Delta y_{ss}}{\Delta u_{ss}}$.
Quickly Simulating a Transfer Function's Response using MATLAB/Octave¶
MATLAB/Octave have built-in numerical integration functions that make simulating a transfer function's behavior very quick. These can all be explored in more detail by typing "help" and then the function name into a MATLAB prompt or an Octave cell. The functions are:
step(), which computes the response of a transfer function system model to a unit (magnitude 1) step inputimpulse(), which computes the response of a transfer function system model to a unit impulse inputramp(), which computes the response of a transfer function system model to a unit ramp inputlsim(), which computes the response of a transfer function system model to an arbitrary input vector.
Many times, using these functions to simulate a system model that you already have in transfer function form is much quicker than writing your own numerical integration solver or solving the system's input-output differential equation analytically (either using the Laplace transform or the method of undetermined coefficients).
Model Construction: Introduction to the Block Diagram¶
- Signals (inputs, state variables, outputs, and combinations thereof)
- Blocks (operators that change one signal into another; usually, these are transfer functions)
- Takeoff points (places where a signal goes more than one place)
- Summing junctions (places where signals are added or subtracted)
We will go through each of these pieces one by one and discuss what each does.
Signals¶
A signal in a block diagram can be an input, an output, a state variable, or some combination of state variables in your system. Usually, a block diagram is drawn with the idea that the signals and operators (blocks) shown are in the Laplace Domain.
Blocks¶
Usually, a block in a block diagram represents a transfer function. If you were to represent the system:
\begin{equation} G(s) = \frac{y(s)}{u(s)} = \frac{1}{s+1} \end{equation}as a block diagram, it would look like this:
Notice that the signals in the diagram are labeled carefully. We say that the block $G(s) = \frac{1}{s+1}$ operates on the signal $u(s)$ to transform it into $y(s)$.
A block is not always a traditional, proper or strictly proper "transfer function" as seen below.
Blocks in 'cascade'¶
When two or more blocks are 'in a row' (or 'in series') in a block diagram, we say that the blocks are 'in cascade.' What this means is that the output from the first block is fed into the second.
A cascade can be replace by one single block that is the product of the two blocks in cascade like this:
Takeoff Points¶
A takeoff point is simply a point where a signal splits and goes to two different places. This concept is illustrated by the graphic below:
Summing junctions¶
At a summing junction, the signal going out of the junction is the sum of the signals coming in.
If it is necessary to subtract one or more signals, a '-' is added next to the corresponding arrowhead.
Blocks in parallel¶
Sometimes, blocks are set up "in parallel" like you might see electrical elements in parallel in a circuit diagram. This requires a combination of a takeoff point and a summing junction.
Disciplined Process: Block Diagram Simplification¶
Why simplify block diagrams?¶
A control system is usually a device or algorithm that is added to an existing system. Because the goal of any control system is to regulate or improve the performance of the plant so it will likely change the dynamics. That is, the controlled system's eigenvalues will likely change. If we do our job right, they will change "for the better."
By simplifying a block diagram we can obtain a new transfer function representing the controlled system and determine its poles and zeros to evaluate whether our controller actually improved its behavior.
Block Diagram Simplification using algebra¶
To simplify the block diagram in the preceding section, the simplest procedure just involves writing equations to describe the relationship between each of the signals. Giving the signals names here is helpful. Write an equation for each of the summing junctions in the diagram, and write an equation that describes all of the blocks in cascade. For our system, we have:
\begin{aligned} e(s) &= r(s)-y(s)H(s) \\ u(s) &= e(s)C(s) \\ y(s) &= u(s)G(s) \end{aligned}By substituting these equations into one another, we get the following:
\begin{equation} (r-yH)CG = y \end{equation}rearranging this equation, we get the following:
\begin{equation} \frac{y(s)}{r(s)} = \frac{CG}{1+CGH} \end{equation}Block Diagram Simplification with multiple inputs¶
Sometimes, it is necessary to look at how a second input to a closed-loop, controlled system might affect the output. For instance, on an autonomous car, a side-wind might affect how well the car tracks the center of the road, even under automatic steering. Therefore, you might want to figure out how bad those effects might be when you design your lane-keeping controller for the car. Often, this type of "disturbance," here labeled $d(s)$, is a common source of a second input to a control system's block diagram. For our example, the side-gust might appear as an additive term on the signal $u$, and might represent the lateral acceleration applied to the vehicle as shown below.
To simplify this diagram, we have to decide which transfer function we're looking for in our final, simplified form.
If we want $\frac{y(s)}{r(s)}$, which shows us how the output follows the "desired" behavior for the system, we would assume that $d(s)=0$ and proceed with the algebraic simplification as above.
If we want to see how our wind gust or "disturbance" $d(s)$ affects the output, we want to look at $\frac{y(s)}{d(s)}$, so we set $r(s)=0$ before simplifying.
Because our system has to be linear for us to be using block diagram algebra and transfer functions to describe it, we can say that superposition applies. That means that for any combination of inputs $r(s)$ and $d(s)$, we can compute the total response of the system using the following simplified block diagram:
To find $\frac{y}{d}$, we have to write the equations for each block and summing junction assuming that $r=0$. This gives us:
\begin{aligned} e(s) &= -y(s)H(s) \\ u_c(s) &= e(s)C(s) \\ u_t(s)(s) &= u_c(s)+d(s) \\ y(s) &= u_t(s)G(s) \end{aligned}Substitution leaves us with:
\begin{equation} (d-yHC)G=y \end{equation}Moving $y$ to one side and dividing appropriately leaves us with the following closed-loop transfer function from the disturbance $d(s)$ to the output $y(s)$:
\begin{equation} \frac{y(s)}{d(s)} = \frac{G}{1+CGH} \end{equation}To find find $\frac{y}{r}$ we would do the same steps, just assuming $d = 0$. Under these conditions the block diagram is identical to the first example so
\begin{equation} \frac{y(s)}{r(s)} = \frac{CG}{1+CGH} \end{equation}note: many example problems of this type are available for you to practice with in Schaum's Outlines of Feedback and Control Systems, and there is also more information in Ch. 6 and 13 of Kulakowski, Gardner, and Shearer's Dynamic Modeling and Control of Engineering Systems.
Graphical interpretation of linear system behavior¶
The keys to understanding a linear system's behavior are:
- The system's eigenvalues, which govern its transient behavior
- The system's zeros, which govern its transient response to an input's derivative(s)
- The system's steady state gain, which governs its steady-state response to an input
Focusing on transient behavior, we learned from the study of linear differential equations that a system's eigenvalues come in two "flavors:"
- A real number $p = -a = \frac{-1}{\tau}$, which contributes a term $Ce^{-at}$ to the system's homogeneous response.
- A complex conjugate pair $p = -\zeta\omega_n \pm \omega_d j$, which contributes a sinusoidal term $C\sin(\omega_dt + \phi)$ to the system's homogeneous response.
By the principle of superposition, all higher-order systems can be broken down into first or second order "pieces," which contribute to the system's overall behavior. This is often accomplished when a system is in transfer function form using a technique like the partial fraction decomposition.
As the number of poles in a system increases, it becomes more and more difficult to understand the system's behavior. To combat this, it is often convenient to plot all of the system's eigenvalues on the complex plane (usually with an 'x' marker), where an eigenvalue's real part is plotted on the x-axis, and its imaginary component on the y-axis. This allows us to quickly ascertain whether the system is stable, whether and how oscillatory it may be, and if necessary, what a particular eigenvalue's quantitative response will be. This is summarized in the plot below.
Interpretation of eigenvalues in the complex plane¶
Complex Numbers¶
In general, the Laplace variable $s$ is a complex number with both real and imaginary components. Up until this point in the course, we have usually only discussed particular values of the number "s" when it is either a pole or a zero of a transfer function. However, there is an entire s-plane out there, and "s" can in theory be anywhere on that plane. We often write a generalized complex number in Cartesian Form as:
\begin{equation} s = \sigma + j\omega \end{equation}Remember that $j$ is defined as $j=\sqrt{-1}$.
We often draw complex numbers on the so-called "complex plane" as vectors, with a real part and an imaginary part.
While imaginary numbers are often plotted like vectors, they are not quite vectors. Adding them and subtracting them might make one think that all rules of vector manipulation apply, but imaginary numbers multiply and divide "normally," whereas vectors do not (they have dot and cross products defined, which are different than algebraic multiplication and/or division). A selection of important properties of complex numbers are listed below:
- To add complex numbers, simply add the real and imaginary components, or use "head to tail" graphical vector addition as shown in the figure above for $s_3=s_1+s_2$.
- To subtract complex numbers, simply subtract the real and imaginary parts of one complex number from those of the other, or use "head to head" graphical vector subtraction as shown in the figure above for $s_4=s_1-s_2$.
- There is no cross product and no dot product for complex numbers. To multiply complex numbers, simply "FOIL" or distribute imaginary and real components. For example, to compute $s_5 = s_1s_2$, one would compute $s_5 = (-3+2j)(1+3j)=(-3)(1)+(-3)(3j)+(2j)(1)+(2j)(3j)=-3-9j+2j-6=-9+7j$.
- It is possible to divide complex numbers, which is another difference between them and 'true' vectors. This is also done algebraically. Often, if the vector is in Cartesian form, it is easiest to multiply the resulting fraction by the denominator's complex conjugate. For example, to compute $s_6 =\frac{s_1}{s_2}$, one might compute: $s_6 = \frac{-3+2j}{1+3j}\frac{1-3j}{1-3j}=\frac{1}{10}(3+11j)$.
Note that MATLAB/Octave can help you compute addition, subtraction, multiplication, and division of complex numbers. MATLAB uses $i$ to denote a complex number rather than $j$ when it displays complex numbers, but it allows you to input $j$ for the complex number if you wish. MATLAB also allows you to extract the real and imaginary components from any complex number. Some examples are shown below.
Alternate forms of complex numbers¶
While it is convenient and familiar to represent complex numbers in Cartesian form as shown above, sometimes it is even more convenient to look at and manipulate complex numbers in either Polar From or Exponential Form.
Polar Form¶
Every complex number has both a magnitude and a direction, just like a vector. To write a vector in polar form, we write:
$$ s_1 = \sigma_1+j\omega_1 = \left|s_1\right|\angle(s_1)$$The magnitude of $s_1$ is computed using the Pythagorean Theorem to find the length of the arrow from the origin of the complex plane to the point $s_1$.
$$\left|s_1\right| = \sqrt{\sigma_1^2 + \omega_1^2}$$The angle of $s_1$ is computed by using the arctangent to find the angle that the line going from the origin of the complex plane to the point $_1$ makes with the real axis. Be careful here-- if your calculator does not have a 4-quadrant "arctan2" function, you will have to use common sense to determine the total angle $s_1$ makes with the real axis. Also note that the angle is not unique, because any addition of $2\pi$ to the complex number's angle would yield the same result.
$$\angle(s_1) = \tan^{-1}\left(\frac{\omega_1}{\sigma_1}\right)$$Exponential Form¶
Writing a complex number in exponential form exploits Euler's formula to represent complex numbers in a very convenient, compact form. Euler's formula states:
$$e^{j\omega} = \cos \theta + j\sin \theta$$Interestingly, this equation gives us both a real and an imaginary component that is derived from "e" raised to an imaginary number. The other thing that you should notice is that this formula represents a unit vector pointing in the direction defined by the angle $\theta$ on the complex plane. This is shown graphically in the figure below:
Feedback Control¶
Anatomy of the closed-loop control system¶
In your last notebook you were introduced to the "canonical" form of a system under feedback control, which is shown in the diagram below:
In the diagram above, as labeled, the system has the following key pieces:
- The plant (labeled G(s) above) is the system, subsystem, or process that is controlled by our feedback controller
- The controlled output (labeled y(s) above) is the output variable of the plant under control of the feedback controller
- The forward path is the transmission path from the sum junction all the way to the controlled output.
- The feedforward elements are all control or compensation elements in the forward path. In the example above, this is just C(s).
- The control signal (labeled u(s) above) is the output signal of the feedforward elements (C(s)) that is applied to the plant G(s).
- The feedback path is the transmission path from the controlled output back to the sum junction.
- The feedback elements (labeled H(s)) relate the actual controlled output y(s) and the feedback signal read by the controller at the sum junction. Most often, H(s) represents a sensor that has some gain and filtering characteristics that apply in the act of measuring the system output.
- The reference input (labeled r(s)) is an external signal signal applied to the system. This is usually our "request" for the system to do something.
- The primary feedback signal (labeled y(s)H(s)) is a function of the controlled output y(s), transformed by the feedback elements H(s). This is the signal that is usually actually read by the controller.
- The error signal is the difference between the reference input and the primary feedback signal.
- Negative feedback implies that the sum junction subtracts r-yH. Positive feedback would imply that the sum junction adds r+yH. Negative feedback is nearly always used in feedback control.
We saw in the last reading that the canonical feedback control system has a closed-loop transfer function:
\begin{equation} \frac{y(s)}{u(s)} = \frac{CG}{1+CGH} \end{equation}But what does this really mean? Note that the controller transfer function, chosen by the engineer, shows up in both the numerator and the denominator of the closed loop transfer function. This means that the addition of the controller and feedback to the plant changes the dynamics of the plant output. This is precisely what we want, if our goal is to improve system performance. In simple terms, $C(s)$ is a transfer function that "makes decisions" about how hard to "push" the plant based on the difference between where the system is at any given instant ($y(s)H(s)$), and where the operator wants the system to be ($r(s)$).
The P(I)(D) Controller¶
In this course, we will focus on variants of the PID or "proportional-integral-derivative" controller, which has the following general transfer function:
\begin{equation} C(s) = \frac{u(s)}{e(s)}= K_p + K_i\frac{1}{s} + K_ds \end{equation}The PID controller responds to the error proportionally (P-control), and/or the error's derivative proportionally (D-control), and/or the error's integral proportionally (I-control). PID control almost always includes a P-control element in the form of $K_p$, but often, $K_i$ and/or $K_d$ can be set to zero depending on the needs of the engineer and/or the plant. This week, we will focus on how variants of this general design work for simple first and second order systems. We'll use simplified examples and the so-called "direct method" of determining closed-loop system performance (see, for more information, chapters 3 and 4 of Schaum's Outline of Feedback and Control Systems or Sections 13.1-13.3 of Dynamic Modeling and Control of Engineering Systems) to build an understanding of how each type of controller works before diving into some more involved design tools that will allow us to apply these controllers effectively to more complex systems.
When synthesizing P(I)(D) controllers, we can have a number of different design goals. Some of these are listed below:
- Above all, we want to ensure system stability, defined as the closed loop system having eigenvalues with real parts strictly less than 0. We also generally want to know how stable a system is.. in other words, how close is it to being unstable?
- If the system is first order, we may want to change its time constant.
- If the system is second order, we may want to specify its natural frequency and/or damping ratio.
- If the system is second order, we may want to specify its rise time, settling time, or percent overshoot.
- We may want the system to achieve a particular performance in tracking a step or ramp input to $r(s)$.
- We may want the system to achieve a particular performance in rejecting an impulse or step disturbance input (not shown above, but disturbances are usually modeled as appearing as an input to a sum junction in between $C(s)$ and the plant $G(s)$).
- We may want the system to achieve particular magnitude and/or phase performance when tracking a sinusoidal input of a known frequency on $r(s)$.
Controller Design: P(I)(D) Controller Disciplined Process¶
- Beginning with a validated system model for the plant, determine the equations and/or system transfer function for each component of the system (e.g. $C(s)$,$G(s)$,$H(s)$, etc.)
- Formulate a model of the total, controlled system by appropriately connecting the system elements using a block diagram. The canonical form above is usually, but not always, a good starting point.
- Determine controller parameters ($K_p$,$K_d$,$K_i$) and determine the system response characteristics under the chosen control paradigm.
Of course, following this design process, you should try to validate your design by implementing it on the real system and see if it performs as expected.
Much of what we'll discuss roughly follows the material in Chapter 10 of Schaum's Outline of Feedback and Control Systems. We will deal with how to handle most of these design goals as needed, but it is helpful to know where we are going. To begin our journey, let's take a closer look at so-called "P-control," which is a special case of PID control with $K_d$ and $K_i$ set to 0.
Controller Design: Proportional (P) Control¶
Proportional control is the simplest form of the PID feedback control paradigm. This controller 'pushes' on the plant proportionally to the error between the reference input and the plant output. As the system gets closer to the desired value of the output, the controller pushes less and less hard. You might think of this controller like a "spring" that drives the plant to the desired output. For a proportional controller, the controller transfer function is a constant $K_p$. Formally, we write:
\begin{equation} \frac{u(s)}{e(s)}=C(s) = K_p \end{equation}Knowing the units of your controller's gain(s) is extremely important. Often $e(s)$ and $u(s)$ will often have different units. This means that the value of $K_p$ will assume the units needed to make a dimensional analysis work out.
For example, if $e(s)$ is a distance error measured in meters, and $u(s)$ is a voltage input to a motor measured in Volts, then the controller's proportional constant $k_p$ must have units of $\frac{V}{m}$.
Proportional-Integral (PI) Control¶
The proportional control term $P(s)$ looks at how far our system is from where we want it to be, but sometimes (for example, for a first-order plant), this is not enough to ensure that our system will be able to track a step input without any steady-state error. The proportional-integral controller solves this problem by looking at how long our system has been away from $r(s)$ as well as how far it is at any given moment. The controller transfer function for a PI controller is:
\begin{equation} C(s) = K_p+K_i\frac{1}{s} \end{equation}This controller, like the proportional controller, gives an input to the system that is proportional to the error signal $e(s)$, but the additional term $\frac{K_i}{s}$ does something different. It pushes on the system proportional to the integral of the error. Think of it this way-- if the system were sitting some constant distance away from the desired value $r(s)$, the integral of the error would increase steadily, and the $K_i$ term would push harder and harder to try to eliminate the error.
Common Controller Design Metrics for Closed Loop Systems Displaying Underdamped Second-Order Behavior¶
While minimizing steady-state error is nearly always a goal of a feedback control system, controlling its transient behavior is also vital to achieving good performance.
Many times, attempting to make a closed-loop system settle quickly to its steady state means that the engineer will need to live with at least some oscillation. Finding the "right" damping ratio $\zeta$ and natural frequency $\omega_n$ for a closed loop system that will act like an underdamped second-order system is often a delicate balance between different design goals.
The two most common design metrics we will use in ME480 for closed loop systems that display underdamped second-order behavior can include, but are not limited to, the following quantities:
Percent Overshoot¶
Often, it is important to limit the system's overshoot in response to a step input. Percent overshoot#:~:text=For%20a%20step%20input%2C%20the,overshoot%20in%20an%20electronics%20context.) can be found approximately for a standard-form, second-order system as:
$$PO = 100e^{-\frac{\pi \zeta}{\sqrt{1-\zeta^2}}}$$This formula can be manipulated to produce a goal for the damping ratio $\zeta$ of the closed-loop control system's dominant eigenvalue pair.
Settling Time¶
Often, the speed at which a system reaches steady state is important. An approximate formula for the $2\%$ settling time of a system that displays underdamped second order behavior can be found as:
$$t_{s,2\%} = \frac{4}{\zeta\omega_n}$$This can be used to help determine a goal for a closed-loop system's damping ratio $\zeta$ if a desired $\omega_n$ is known, or to determine $\omega_n$ if the goal for $\zeta$ is known.
Proportional-Derivative (PD) Control¶
The transfer function for PD control is, in general,
\begin{equation} C(s) = K_p + K_ds \end{equation}In contrast to the PI controller, which reacts to both the magnitude of the error the area under the error curve $e(t)$ as the system attempts to reach steady state, the PD controller reacts to the error through $K_p$ and the speed at which the error signal is changing through $K_d$, which reacts proportionally to the error signal's time derivative. You might think of this as a type of "virtual damper" the way you might think of $K_p$ as a virtual spring. PD control is very popular for second and higher order systems. At this point in the course, we're still limited to analyzing the behavior of a closed-loop system by directly studying the solutions to its characteristic equation, so let's look at a second order example.
Controller Design: Root Locus Method¶
Up until this point in the course, we have been working with the following "canonical form" of a feedback control system:
To choose controller gains for our controller $C(s)$, we have been directly manipulating the control system's closed loop characteristic equation after finding it using the system's closed loop transfer function:
$$\frac{y(s)}{r(s)} = G_{cl}(s)=\frac{CP}{1+CPH}$$Again, when we haved designed a controller to particular performance specification, like $\zeta$, we are implicitly setting a goal for the system's closed loop eigenvalues. Because manipulating third and higher order polynomials is difficult, it is hard to determine the closed loop eigenvalues that will yield the desired closed-loop behavior. So we will turn to the new Root Locus Plot tool. This tool will be able to reveal how to choose our controller's structure (P, PI, PD, or PID) and set its gains to achieve a desired system eigenvalue.
The Root Locus Plot is a graphical description of how a system's closed loop eigenvalues move around the s-plane when one controller parameter (e.g. $k_p$,$k_i$,etc) varies from $0\rightarrow\infty$. This might seem restrictive, because it won't immediately allow us to choose all gains in a PID controller, but we have to crawl before we can walk! Once we learn rules for quickly sketching root locus plots, we'll be able to use them to determine:
- What eigenvalues are and are not possible with a particular controller, particular plant P(s), and sensor transfer function H(s).
- How to set controller gains to achieve eigenvalues that are possible with our controller
- A measure of "how stable" our control system is for a given controller design.
Preparing your system for root locus analysis¶
To begin, we will assume that we can factor one parameter out of our controller transfer function $C(s)$. After factoring this one, real "gain factor" out of $C(s)$, we will denote the "static" part of our controller transfer function $C^\star(s)$. This is shown in the figure below:
Determining the Closed Loop Characteristic Equation¶
In the closed loop transfer function shown above,
$$\frac{y(s)}{r(s)}=G_{cl}(s) = \frac{KG}{1+KGH}$$,
The characteristic equation for the closed loop system is found by looking at the denominator of $G_{cl}(s)$. In other words, the characteristic equation is given by:
$$1+KGH=0$$.
Instead of "clearing fractions" in this characteristic equation, the way we did when we were manipulating it directly to achieve particular eigenvalues, we will leave the transfer function $GH$ as a fraction of polynomials. In some texts, like Schaum's outline, the authors explicitly say $G(s)H(s) = \frac{N(s)}{D(s)}$ to denote that it is a (proper or strictly proper) fraction of the two polynomials $N(s)$ and $D(s)$. We won't worry too much about that. What we will do, however, is rearrange the characteristic equation using some simple algebra to obtain:
$$-\frac{1}{K}=G(s)H(s)$$.
This is important because it means that because $K$ is a positive real number, the so-called "open loop transfer function" (this is an unfortunate name, but that's what people call it... see the appendix) has to reduce to a negative real number for any value of $s=\sigma+j\omega$ if it is a solution to the system's closed-loop characteristic equation.
That's it. That's the whole key to learning to draw, interpret, and use root locus plots in control system design. Because $-\frac{1}{K}$ is a negative real number, $G(s)H(s)$ must be a negative real number also when evaluated at a point $s_1$ on the complex plane that is a solution to the characteristic equation.
Finding Eigenvalues on the Root Locus and their Gains¶
Let's look a little closer at what the realization above tells us. If we were to plot the point $s=-\frac{1}{K}$ on the complex s-plane, we would get something like this:
The angle of this point will always be $\pm180^\circ$, possibly with $360n^\circ$ added without loss of generality. This leads to what we call the "angle criterion" for the root locus plot.
Angle Criterion: Determining Eigenvalues on the Root Locus¶
A point $s1=\sigma_1+j\omega_1$ is "on the root locus" of a system if the angle of the transfer function $G(s)H(s)$ is $\pm 180^\circ$ when that transfer function is evaluated as a complex number at $s=s1$. The angle criterion can be written formally as:
$$\angle-\frac{1}{K}=\angle G(s)H(s) = \pm 180^\circ$$The angle criterion is what we'll use to derive all of the rules for sketching root locus plots by hand. It does not tell us how to get to that point by setting the gain factor $K$ to a particular value. It merely tells us that we can get to that point, which means that if the point $s_1$ is a desired eigenvalue for our closed-loop system, we'll be able to know whether it's possible to achieve it, regardless of how complex the transfer function $G(s)H(s)$ is! If you need help figuring out what "the angle of G(s)H(s)" means, please see the other notebook in this folder, which is a review of complex numbers.
The "magnitude criterion" does tell us how to get to a particular point on the root locus. As you might guess, if the closed-loop characteristic equation, which we write as $-\frac{1}{K}=G(s)H(s)$, is satisfied, then we can assert that the magnitude of $-\frac{1}{K}$ must equal the magnitude of the transfer function $G(s)H(s)$.
Magnitude Criterion: Determining Gains for Possible Eigenvalues¶
If a point on the root locus is given by $s=\sigma+j\omega$, then the following equality holds:
$$\left|-\frac{1}{K}\right|=\left| G(s)H(s)\right|$$Therefore, we can find the gain factor $K$ that 'gets us to' a particular point $s_1$ on the root locus for our system using $K=\frac{1}{\left|G(s)H(s)\right|}$.
What is the meaning of this "open loop transfer function" G(s)H(s)?¶
One of the most common points of confusion for students learning to draw the root locus is this idea that we call $G(s)H(s)$ the "open loop transfer function." It's not a very good name, but it is convention. In reality, the "open loop transfer function" as it's defined in root locus analysis is a very specific transfer function indeed:
<a href=" ../W07A_Control_ComplexPlane_RootLocusFundamentals/W07A_Control_ComplexPlane_RootLocusFundamentals.html#What-is-the-meaning-of-this-"open-loop-transfer-function"-G(s)H(s)?" > Link to Original Context in: W07A_Control_ComplexPlane_RootLocusFundamentals</a>
Summary of Root Locus Properties¶
- The characteristic equation for this system is $1+KG(s)H(s)=0$, which can be rewritten as $-\frac{1}{K}=G(s)H(s)$.
- The location of each possible pole is governed by the angle criterion for the root locus: $\angle -\frac{1}{K}=\angle G(s)H(s) = \pm 180^\circ$.
- There are exactly as many branches (paths for a closed loop pole to follow) as the number of poles in G(s)H(s).
- When $K=0$, the closed loop poles lie at the poles of the 'open loop transfer function' G(s)H(s) because of the magnitude criterion: $\left|-\frac{1}{K}\right| = \infty = \left|G(s)H(s)\right|$. The only way to make $G(s)H(s)=\infty$ is to set $s$ to one of the poles of $G(s)H(s)$.
- When $K=\infty$, the closed loop poles go to either a zero of $G(s)H(s)$ or infinity. This is because of the magnitude criterion. The fraction $\left|-\frac{1}{K}\right|=0= \left| G(s)H(s)\right|$. There are two ways to make $G(s)H(s)=0$. One is that $s$ is one of the zeros of G(s)H(s), and the other is that the magnitude of $s$ goes to infinity. Assuming that the transfer function is strictly proper (meaning that at least one branch of the root locus goes to infinity), setting the magnitude of $s$ to infinity will make the transfer function's magnitude go to zero.
- If a branch goes to infinity, it does so along an asymptote.
- A branch will lie on a portion of the real axis if there is an odd number of poles and zeros of $G(s)H(s)$ to the right of that portion. This is required by the angle criterion.
- The root locus branches can "break away" from and "break in" to the real axis. They always do so perpendicularly to the real axis. This is required by the angle criterion.
- The root locus is always symmetrical about the real axis. This is required by the angle criterion.
Disciplined Process: Steps 1-5 of sketching the root locus¶
- Find the Closed loop characteristic equation for your system: $$1+KG(s)H(s)=0$$
- Write the Closed loop characteristic equation as $$\frac{-1}{K}=G(s)H(s)$$
- Plot the poles and zeros of $G(s)H(s)$ on the s-plane. Identify any segments of the R.L. on the real axis: the root locus lies on a portion of the root locus if and only if an odd number of "open loop" poles and zeros (poles and zeros of $G(s)H(s)$ are to the right of that portion of the real axis.
- Identify the angle(s) of any asymptotes and the point on the real axis where the asymptotes intersect: $$\sigma_a = \frac{\sum poles -\sum zeros}{n-m}$$ $$\theta_a = \frac{2k+1}{n-m}\pi$$
- If you suspect that the root locus crosses into the right half of the s-plane (based on the asymptote angles), calculate the crossing frequency (where the system crosses into unstable territory) and critical gain (gain that makes this happen) using the closed loop characteristic equation. Plug in the point $s=s_{crossing}=j\omega_{crossing}$ into the open loop transfer function $G(s)H(s)$ to find the value of $K=K_{critical}$ that will result in an unstable system.
Finally, use the asymptotes to sketch each branch (one per pole of $G(s)H(s)$) of your system's root locus, noting that each pole of the closed loop system will begin at one of the poles of $G(s)H(s)$ and terminate at either a zero of $G(s)H(s)$ or at infinity. Approximate any break-in or break-out points using the real axis rule (step 3).
Controller Design: The Root Locus Angle Deficiency Method¶
At this point, you probably have a pretty good handle on what it means to sketch a root locus, but you might be wondering how it is used to design controllers for systems that may be too complex for the "direct method" of manipulating the characteristic equation.
The process we'll follow is actually very similar to the process we used for calculating angles of departure and arrival in our root locus sketching steps. This design process is called the angle deficiency method of controller design. The design process consists of a few simple steps:
- Sketch the root locus of our system under proportional control, along with a target "dominant" eigenvalue $s_d$ that we wish the system to have after the design is complete. If the rightmost (slowest) branches pass through a suitable eigenvalue that meets your design requirement, use the magnitude criterion to find the necessary proportional gain $K_p$ to achieve your goal, and proceed to step 6. If not, proceed to step 2.
- If the system is unable to reach our performance target under proportional control, determine whether adding a first-order zero (PD control), a first-order zero and a pole at the origin (PI control), or a second-order zero and a pole at the origin (PID control) might "bend" the root locus in the proper direction to allow you to hit your desired eigenvalue location. In general, if you need to "bend" the rightmost root locus branches to the left, try PD control first, and if you need to "bend" the rightmost branches to the right, try PI control first.
- Using PD, PI, or PID control, and the angle criterion for the root locus, calculate the location of the controller zero(s) that will allow the root locus to hit your desired eigenvalue location $s_d$.
- Using the magnitude criterion, calculate the "root locus gain" $K$ that will allow you to achieve your desired eigenvalue with your chosen zero location.
- Using your controller transfer function, the now-known zero location(s) and root locus gain, calculate the unknown controller gains needed to implement your controller.
- Simulate the response of your system under closed-loop control given your new controller design. Pay particular attention to the minimum and maximum controller outputs (plant inputs) that result from a step input on $r(s)$ of a magnitude relevant to your overall design goals. If the controller is physically feasible, you're done! If it is not, you may have to rethink your design.