Newsletters - MATLAB Digest
How to Build a Clock or Controlling an Oscillation in a Nonlinear System Using MATLAB®, Simulink®, and the Control System Toolbox
by Carla Schwartz and Richard Gran, The MathWorks, Inc.
I. Abstract
The MathWorks, Inc. tools include MATLAB, Simulink, and the Control System Toolbox. These tools provides a substantial collection of analysis and design methods for control systems. The Control System Toolbox has been improved; some of the latest features of Versions 4.0, 4.1, and 4.2 of the Control System Toolbox include:- Frequency Response Data (FRD) models that use frequency responses to model systems
- The ability to collect a set of LTI models into an array under a single variable name
- The LTI Viewer for displaying time and frequency responses of linear systems
You can use Simulink to simulate the behavior of both linear and nonlinear models, and an explicit LTI object may be included in the model.
In this article we use these tools to perform frequency domain analysis on a nonlinear model of a mechanical pendulum clock. A Simulink model that creates a describing function for an arbitrary nonlinear system is developed and is applied to the design of the clock escapement. We use the describing function of the escapement to find the angle for the escapement gear (interpreted as a a control gain that translates into a slope for the escapement gear) that causes the oscillation of the pendulum to fit within the clock case. The escapement is designed to provide sufficient energy to the pendulum so that wear on the bearings does not cause the clock to stop.
The analysis proceeds in four steps:- We use the describing function of the escapement to develop an array of FRD models, and, in the process, develop a Simulink model to create the arbitrary describing function of an arbitrary nonlinearity.
- We analyze the dynamics of the pendulum and build an array of FRD models for the pendulum that can be used with the array of escapement FRD models.
- We use the LTI viewer to analyze the frequency response of the array of FRD models. From this analysis, we calculate the gain that causes the system to oscillate with a specified amplitude for the pendulum's swing.
- We simulate the design using Simulink to validate the analysis results.
The files used in this paper are available in the clock.zip file.
II. Introduction
The mechanism for a pendulum clock is shown in Figure 11.
![]() |
Figure 1. A Clock Pendulum Mechanism, Including the Pendulum (1 to 8), the Escapement Wheel (9), and the Anchor (13) 1 This copyrighted image is reprinted with the permission of the British Horological Society. |
- It provides energy to the pendulum so the clock doesn't stop as a result of friction.
- It regulates the speed of the gear motions to the natural period of the pendulum so that the timing of the clock is accurate.
- It provides the feedback mechanism needed to ensure that the speed of the clock is not affected by the force from the driver (springs or weights) and wear in the various mechanical parts does not affect the operation of the clock.
This type of mechanism was responsible for the emergence of accurate mechanical clocks in the eighteenth century. The mechanism relies on a mechanical invention called an "escapement" (the parts labeled 9 through 15 in the diagram). The escapement shown was invented by Christiaan Huygens in 1656. With this mechanism, he was able to build a clock that was accurate to within 10 seconds per day. By modern standards this may not seem remarkable, but previous clock mechanisms did not have the regulation capability of the Huygens escapement. They would slow down as the springs unwound or as the weights that drove the clock dropped down. Earlier clocks could vary by as much as several hours if they were not reset daily.
Before we go into the details of the design of the escapement, let's look at how it works. The pendulum is mounted on a rigid structure in the clock that is called a "back cock" (Figure 1, part 1). This allows the pendulum to swing freely, unencumbered by the mechanism that drives the hands and the escapement. The pendulum suspension (part 2) is carefully designed to allow the pendulum to swing with minimal friction. The escapement imparts its force to the pendulum through the crutch (part 3) and the crutch pin (part 4). The pendulum itself (parts 5 and 7) moves back and forth with a period that is determined solely by its length (as was shown by Galileo Galilei in 1582). The pendulum includes a nut (part 8) at its bottom that can be adjusted to change the effective length and thereby adjust the period of the pendulum.
The escapement wheel (part 9) is linked to all of the gears in the clock (both the hands and the windings) through gears that are not shown in the figure. The escapement wheel is what we will be designing in this paper.
For the remainder of this description, we will focus on the escapement anchor (part 13) and its two sides. On the left side of the anchor is a triangular face called the entrance pallet (part 12) and on the right side is the exit pallet (part 15). As the pendulum moves to the left, it forces the anchor to rotate and ultimately disengage its face (part 11) from the escapement wheel (part 9). When it does, the exit pallet (part 15) engages the next tooth of the escapement wheel (part 9) and the gears move. If the pendulum is set to swing at a period of 1 second, the result of the motion of the escapement gear is to move the hands on the clock by 1 second. When the right face of the anchor engages the gear at the exit pallet (part 15), the anchor (part 13) and the escapement wheel (part9) interact. The escapement wheel (part9), which is torqued by the winding mechanism, applies a force back through the anchor and then through the crutch and crutch pin (parts 3 and 4) to the pendulum. The alternating release and reengagement of the anchor and escapement wheel provides the energy to the pendulum that compensates for the energy lost to friction. The distinctive and comforting "tic-toc" sound that a clock makes is a by-product of this interplay of the pendulum and the escapement wheel.
In this article, we use MATLAB, Simulink, and the Control System Toolbox to design the teeth on the escapement wheel of a clock. As a consequence of this design, an escapement not only helps to regulate the clock but also compensates for wear in the mechanism of the clock.
III. A Simulink Model for the Clock
The clock model has two parts: the model of the pendulum and the model of the escapement. The pendulum model comes from the application of Newton's laws by balancing the force due to the acceleration with the forces created by gravity. If the pendulum angle is![]()
Equivalently, substituting for J, the model is:
![]()
These dynamics are simulated in the pendulum clock Simulink model, Clocksim, shown in Figure 2. This model was built using the Simulink primitive blocks for integrators, gains, a summing junction, and the nonlinear sine block.
|
Figure 2. The Simulink
Model for the Pendulum Clock
Click on the graphic for the full size of the image. |
The angular acceleration of the pendulum, Theta2dot, is formed at the right of the sum block, and is calculated using the following three terms:
- The damping acceleration from the equation above (dampcoef * Thetadot)
- The gravity acceleration from the equation above (g * sin(Theta))
- The acceleration created by the escapement
The coefficients (J, m, l, c (dampcoef) and g) are introduced through the M-file clockdat.m when Clocksim is loaded.
The Simulink model creates the variables representing the angular position and velocity, Theta and Thetadot, by integrating the angular acceleration Theta2dot. In addition, an initial pendulum angle of 0.15 radians (about 8.6 degrees) is provided to start the clock (a pendulum clock only sustains its oscillation if the pendulum is moved to an angle that causes the escapement to operate).
A Simulink model for the nonlinear escapement mechanism is shown in Figure 3 below.
|
Figure 3. Simulink Model
of the Escapement Dynamics
Click on the graphic for the full size of the image. |
The attributes of the escapement that are modeled here include:
- The angle of the pendulum at which the escapement engages (applies an acceleration) is denoted by Penlow (as can be seen it is given a value of 0.1 radians)
- The angle of the pendulum at which the escapement disengages (jumps from one gear tooth to the next) is denoted by Penhigh (with a value of 0.15 radians)
- The force is applied only when the pendulum is accelerating (i.e. when the pendulum is moving toward the center of its swing). You know the pendulum is accelerating when the derivative of the pendulum angle is positive.
- Since the pendulum motion is symmetric, the calculations are only made for positive angles (see the absolute value function block in Figure 3). The force applied is given the correct sign by multiplying it by +1 or 0, depending on the sign of the angle.
Design Parameter for the Escapement Model
The controller design gain for this problem represents the level of
acceleration that needs to be applied when the escapement is torquing
the pendulum. This ultimately depends on the design of the escapement
teeth (parts 10 and 6 in Figure 1). This acceleration is determined
from the force applied by the escapement, and it appears in the Simulink
model of Figure 2 as the gain, Design Gain. The units of this gain
are rad/sec2/unit input, or length over inertia. Another
attribute of the escapement, also modeled in Figure 2, is the spring
dynamics associated with the escapement. These dynamics can be quite
complex in their details. Fortunately the clock doesn't require that
these dynamics be modeled accurately. A simple first-order time constant
(0.01 seconds) is used to model the escapement spring dynamics.
IV. The Escapement Control Design Problem Statement
The design problem we are trying to solve here is to choose the value of Design Gain that causes the pendulum system (the nonlinear model above) to oscillate. This control gain must be adjusted so that the system is unstable (i.e., oscillates). The amplitude of the oscillation must also be small enough so that the pendulum doesn't hit the walls of the clock. We select an oscillation amplitude of 0.24 radians (or about 14 degrees) as the desired value. This value is arbitrary, but it must be greater than the value Penhigh in order to have the escapement work. The gain must also be sufficiently large so that the pendulum sustains the oscillation while the driving force drops as the spring on the clock unwinds.
V. Describing Function Analysis
The describing function for a nonlinear system is a linear gain that depends on input amplitude and frequency.
The function is computed by assuming that the input to the nonlinearity
is a sine with amplitude E and frequency
.
The output of the nonlinearity is then a periodic function that can
be analyzed using a Fourier series. The first sinusoidal component
(the fundamental) of this series is used to define the gain. It is
the result of dividing the fundamental by the input amplitude.
Describing function methods have traditionally been used as a first-order modeling technique for predicting oscillations in systems with nonlinearities. See [2], [3], and [4] for references on describing functions. The describing function is the amplitude (and sometimes frequency) dependent gain obtained using the first harmonic of the response of the nonlinearity to a set of sinusoids at various frequencies and amplitudes.
The describing function depends on the input amplitude and frequency, but is otherwise simply a gain, so hence, it can be treated as a linear gain for the analysis. This gain provides an approximation to the behavior of the nonlinear system. With it, you can analyze the first-order behavior of the system's response to a sinusoidal input. Since the gain has an implicit value determined by the amplitude of the input, when oscillatory behavior of the system is determined, the amplitude of the oscillation will also be determined. Our goal is to find the clock input amplitude (angular position) that predicts a sustained oscillation at a desired amplitude.
When a nonlinearity is embedded in a system, it may seem odd that sinusoidal inputs are assumed. This assumption is justified when the input to the nonlinearity is filtered in such a way so that the higher frequencies at the output of the nonlinearity are negligible. In the case of the clock model, this is a valid assumption due to the remaining pendulum and escapement spring dynamics. This makes the clock design problem ideal for the application of describing function techniques.
To sustain a stable oscillation in a nonlinear system, the describing function gain must be such that the open-loop transfer function has a maximum gain of 0 dB and a phase shift of 180 degrees at the frequency of oscillation. We use this as the condition when we analyze the clock and predict the pendulum's frequency of oscillation and amplitude.
VI. Frequency Response Data Models in the Control System Toolbox
We've posed the problem of controlling the oscillation in a pendulum clock. We want to adjust the value of the variable Design Gain in Clocksim so that the amplitude of the oscillation is large enough to sustain the oscillation, but not so large as to exceed the physical limitations of the size of the clock. For this example, we want an oscillation amplitude of about 0.24 radians, or 14 degrees.
The clock design problem is perfect for the application of the frequency domain analysis functions of the Control System Toolbox. These functions operate on special structures in the toolbox called Linear Time Invariant (LTI) objects. You can use LTI objects to store data for LTI models as transfer function models, zero-pole-gain models, or state-space models. For example, you can apply the bode command to an LTI object called sys, to display the Bode plot for the LTI model sys.
There is a new type of LTI object for modeling linear systems in Version 4.2 of the Control System Toolbox. This is the frequency response data (FRD) objectmodel. An FRD object model stores frequency response data as a table that consists of the complex amplitude at specified frequencies. If you only have frequency response data available, you can use FRD models to represent linear system behavior in a feedback control system. As is the case for all LTI objects, the FRD object is overloaded so you can apply the usual operations such as *, /, +, and - to FRD models.have the meanings that one would expect.
Since describing functions are determined by the amplitude-dependent response data of the nonlinear system obtained from a Fourier series decomposition. Thus, the describing function for a given nonlinear system they can be used to represents an FRD model for a given nonlinear system.
VII. LTI Arrays of Models
Version 4.2 of the Control System Toolbox has a new feature called LTI arrays. You can use this feature to collect a family of models in an array under one variable name.
Most of the Control System Toolbox operations (such as addition and multiplication) and many of its functions (such as feedback and series) are overloaded for LTI arrays. This means that you can add two LTI arrays containing the same number of models, and each model in the resulting LTI array is the result of the sum of the corresponding models from the original LTI arrays.
For each amplitude and frequency of the sinusoid that we apply to the nonlinear system modeling the pendulum clock escapement, we obtain a set of frequency response data. We can collect these amplitude-dependent sets of frequency response data in an LTI array of FRD models. We do this in section IX.
First, let's compute the describing functions and collect the necessary frequency response data for constructing an array of FRD models to represent the escapement nonlinearity.
VIII. Computing Describing Functions for Nonlinear Models
Mathematically, the describing function of a nonlinearity f(x) is defined as follows:
Assume the input to a nonlinearity f(x) is a sine with amplitude E and
frequency
(i.e. a period
of T = 2
/
).
![]()
The describing function DF is computed as

You can use a combination of MATLAB and Simulink to compute the describing function of a nonlinear model. The model, GenerateDescribingFunction, is shown in Figure 4. This model implements the mathematical definition of the describing function presented above using a vector of amplitudes, E, that are in the gain block titled "Amplitude Vector from Workspace." The model in Figure 4 calculates the describing function for all of the defined input amplitudes and for a single frequency. Vector calculations are denoted using the wide (vector) lines in the model. This model is executed from within the function GetDF.m.
|
Figure 4. Simulink Model
Used to Calculate Describing Functions
Click on the graphic for the full size of the image. |
You can use the model in Figure 4 and code to generate the describing function for any nonlinearity you choose. To do this:
- Substitute your model for the block labeled Nonlinearity Under Test in the model GenerateDescribingFunction.
- Specify the vector of amplitudes, E and frequencies, omega in the file GetDF.m.
- Type GetDF at the command line.
For the purpose of the clock analysis, the block labeled Nonlinearity Under Test is the same as the nonlinear escapement model depicted in Figure 3.
Notice that the frequency response data generated by GetDF.m is stored in a multidimensional array called array. The dimensions of array are 1-by-1-by-length(omega)-by-length(E). The first two dimensions represent the number of outputs and inputs. Each entry in array gives the describing function gain (complex) for a given frequency (in omega) and input amplitude (in E). This is precisely the format needed to store the frequency response data in order to create an array of FRD models. Notice that the dimension corresponding to frequency immediately follows the dimensions corresponding to the outputs and inputs. In this example, the length(E) is 29. This results in 29 FRD models in the array that we construct in section IX, one model corresponding to each input amplitude.
The result of running GetDF.m for the escapement nonlinearity is shown in Figure 5. The amplitude and phase are plotted using the code in GetDF.m.
|
Figure 5. Describing Function
Amplitude (Left) and Phase (Right), Plotted As a Function of
the Amplitude of the Input and the Frequency of the Input Click
on the graphic for the full size of the image.
|
This particular nonlinearity does not have a fundamental Fourier component that varies with frequency, so the describing function gain is constant for all frequencies. It does, however, have a phase shift. The describing function for this nonlinearity was determined analytically in [2], and the values we calculate here are precisely as those derived there.
Describing Functions and Oscillations
The describing function gain given in the magnitude plot above starts
at zero when the amplitude is zero, and increases to a maximum value,
after which it decreases. For clarity, only the value of the gain
after the maximum occurs is plotted in Figure 5. This is the region
of operation of the clock.
Let's examine why the system operates in this region of the gain plot. When the system gain is set so that the system oscillates, the amplitude of the oscillation is at some value (say Eo). If the system were to change so that the effective gain decreased (for example, if the spring force changes), then in this region of the frequency response curve, the describing function gain would increase. This is because a decrease in amplitude causes an increase in gain (Figure 5). In a similar way, if the system gain were to increase (e.g., after winding the clock), then the amplitude increases, and the describing function gain decreases. This reasoning implies that not only is the oscillation sustained but that the oscillation is stable (in the sense that it does not increase in amplitude) as well. For similar reasons, the values of the describing function gain before the peak in the curve are not in the stable region of oscillation of the pendulum.
IX. Frequency Response Analysis with the Control System Toolbox
A. Creating the Array of FRD ModelsWe have already stored the frequency response data in the format required for constructing an array of FRD models in the variable array, when we ran GetDF.m. To construct an array of FRD models Dfsys from the multidimensional array of complex frequency response data stored in the variable array, along with the frequency response vector omega, type
Dfsys = frd(array,omega);
This creates an LTI object that is an array of FRD models. The FRD represents the amplitude-dependent set of linearized frequency response models for the nonlinear dynamics of the escapement shown in Figure 3.
To analyze the complete system, we must include the remaining linear
pendulum and escapement spring dynamics shown in Figure 2 (using the
approximation sin(u)
u).
We create a single-input, single-output (SISO) transfer function model
of the pendulum called dynamics. We then multiply each model in the
FRD array, Dfsys, by the transfer function for the model dynamics.
The product is a new FRD array called clocksys, and it is the transfer
function of the open-loop pendulum clock model.
A block diagram for this model is shown below in Figure 6 and the M-file is FRD_Analysis.m
|
Figure 6. Block Diagram
of the Composite Models for the LTI Array of FRD Models "clocksys" in
the FRD_Analysis.m file. Click to enlarge image.
|
This M-file also includes the command for opening the LTI Viewer with the Bode plot of the FRD array clocksys.
B. Analysis
The LTI Viewer is a graphical user interface (GUI) that provides a
broad set of linear analysis plotting functions. When an array of
LTI models is imported into the LTI Viewer, you can see the plots
of the entire collection of models displayed simultaneously. You
can only use the frequency response plotting functions such as bode
or nyquist to analyze FRD models.
We begin the analysis of the (linearized) frequency response of the clock model by examining the Bode plots of all the models constructed from the describing function frequency response data. We are trying to predict the amplitude of oscillation of the clock and eventually ensure its arc of oscillation is about 0.24 radians.
Recall that a linear system maintains a stable oscillation (neither increasing or decreasing in amplitude) at a particular frequency if the maximum gain of the Bode plot at that frequency is 0 dB, while the phase at that frequency is ± 180 degrees. This is because when a linear system has a transfer function that is 0, the output is fed back exactly in phase with the input. The nonlinear escapement adjusts itself so that the loop gain is maintained at exactly -1 through changes in the oscillation amplitude. Our first objective is to see if we can predict the oscillation amplitude from the collection of Bode plots derived from the describing function data.
The plot shown in Figure 7 comes from loading the Bode plots of the array of FRD models clocksys into the LTI Viewer using the command:
ltiview('bode',clocksys)
|
Figure 7. The LTI Viewer
Displays the Bode Plots of All 29 Models in the LTI Array clocksys.
Click to enlarge image.
|
X. The Right-Click Menu in the LTI Viewer
Another new feature of the LTI Viewer is the ability to control many features of the plots with a menu that is invoked using the right mouse button (right-click menu). This menu includes an option for opening a GUI and selecting a subset of models in an LTI array whose plots you want to display. To access this menu, just right-click anywhere in the plot region of the LTI viewer away from the actual plot lines. Figure 8 shows how to open the LTI array Model Selector GUI. Figure 9 shows the resulting GUI.
|
Figure 8. Opening the
LTI Array Model Selector from the Right-Click
Menu of the LTI Viewer. Click to enlarge image. |
|
Figure 9. The LTI Array
Model Selector. Click to enlarge image.
|
XI. Selecting Models in the LTI Viewer Using the LTI Array Model Selector GUI
We can select models with the Model Selector GUI using one of two methods: •Indexing into the dimensions (selecting entries in the first column of the array of models) •Selecting models based on bounds you assign to the plot characteristics
Since we want to see if we can select only that model whose peak magnitude response is 0 dB, while the phase is ± 180 degrees, we select models using the Bound on Characteristics method. This selection requires a special syntax that is shown in Figure 10. In this example, we select the only model whose peak gain is between ± 1 dB.
|
Figure 10. Selecting Models
Using the Bound on Characteristics Method.
Click to enlarge image. |
Once you press the Apply button, the LTI Viewer only displays one plot. You can open the right-click menu once more and select Peak Response under the Characteristics submenu (Figure 8) to mark the peak response of the magnitude plot. Then right-click directly on the peak response marker and you can see the plot data for that point, as shown in Figure 11.
XII. Checking the Simulation Results Without Performing Gain Design
Figure 11 shows the Bode plot of the 19th model in the array of FRD models. The frequency of oscillation is determined by the frequency at which the phase angle is 180 degrees. The combined effect of low damping in the pendulum and the escapement spring dynamics ensures that this frequency is very close to the pendulum natural frequency.|
Figure 11. The LTI Viewer
Displaying the Bode Plot of One Model and Its Peak Response
Click to enlarge image. |
Note that the frequency of oscillation is very close to the resonant frequency of the pendulum because of the rapid phase shift. The phase goes from about 20 degrees to almost 200 degrees at the resonant frequency because of the light damping in the pendulum.
We don't quite have 0 dB for the peak response of the gain for this model, since we didn't supply values of the input in the array E that exactly matches the 0 dB condition. By checking the value of E(19) in MATLAB, we see that this model corresponds to an input amplitude of 0.65 radians. The model corresponding to E(20) in the array has a peak response slightly less than 0 dB. It represents an oscillation at an amplitude of 0.7 radians. At this value, the gain is too low to sustain the oscillation, causing the amplitude to decrease. This decrease increases the gain (modeled by the describing function), and causes the system to again oscillate. In other words, the operating point is stable.
From this analysis, we expect that the amplitude of oscillation of the pendulum would be somewhere in this vicinity. Let's simulate the actual nonlinear clock model (with Design Gain set to 1) and see what the amplitude of oscillation is.
The (zoomed in) results of this simulation are shown in Figure 12. Notice that the amplitude of oscillation is where we expected it to be from the analysis.
|
Figure 12. Zooming in
on the Results of Simulating Clocksim Without Designing the
Gain Amplitude of Oscillation (Pend. Angle) Is Approximately
0.7 Radians.
Click to enlarge image. |
XIII. Selecting Design Gain in the Simulink Model For a 0.24 Rad. Oscillation
So far we've shown that for a Design Gain of 1, the clock model, Clocksim, shown in Figure 2 predicts that the pendulum oscillates at about 0.7 radians.
This oscillation is too large, so we need to reduce the force applied to the pendulum. We need to change Design Gain so that the oscillation occurs at 0.24 radians. You might think this is a simple matter of scaling. However, we need to take the nonlinearity of the escapement model (Figure 3) into account. To do this, let's return to the LTI Viewer and examine the response of the model corresponding to an input amplitude of 0.24 radians. This input amplitude provides too high a gain so it does not correspond to a stable sustained oscillation.
First note that the model that corresponds to an oscillation amplitude of 0.24 radians is the 10th model in the array clocksys. You can see this by typing E (the name for the amplitude variable) at the command line.
Now we can use the indexing into the dimensions method for displaying the Bode plot of the 10th model in the LTI array of FRD models. This model selection method and its resulting display are shown in the Figures 13 and 14.
|
Figure 13. Selecting the
10th Model by Typing in the Textbox.
Click to enlarge image |
|
Figure 14. The Bode Plot
of the 10th Model in clocksys, with Its Peak Response Displayed
Click to enlarge image |
In this figure, we used the Characteristics submenu from the right-click menu to show the peak response, just as we did in Figure 11.
Notice that the peak amplitude is 18.8 dB. If we scaled this model, we could force the peak response to be 0 dB. The required scaling factor is:

This is what we set the Design Gain to in Figure 2.
XIV. Simulation Results Using the Design Gain
After changing Design Gain in Clocksim and rerunning the simulation (Figure 15), we see that the clock does indeed oscillate with an approximate arc length of 0.24 radians, and we have designed the escapement to our specifications.
|
Caption goes here.
|
XV. Conclusions
We have presented a Simulink model for a pendulum clock, along with a Simulink model for describing function analysis. We used the latest features of the Control System Toolbox to:
- Create the describing function
- Create an amplitude-dependent array of frequency response data (FRD) models
- Analyze the Bode plots in the LTI Viewer
- Analyze the simulation results from the Simulink model
- Design the escapement gain so as to achieve a desired stable oscillation amplitude for the pendulum
Our analysis using the LTI Viewer included:
- Selecting models based on their peak response
- Displaying peak response characteristics on the Bode plot
- Determining the correct Design Gain for the desired amplitude of oscillation
- Selecting models based on their index in the array of models
Bibliography
[1] British Horological Institute, "Know your terminology - clock
time train," [http://www.bhi.co.uk/hints/clock.htm],
accessed November 19, 1998.
[2] Gibson, J.E., Nonlinear Automatic Control, McGraw-Hill,
1963.
[3] Slotine, J-J.E., and W. Li, Applied Nonlinear Control,
Prentice Hall, 1991.
[4] Rimer, M. and R. Gran, "Analysis of Systems with Multiple
Nonlinearities," IEEE Transactions on Automatic Control,
January 1965.
Store















