The main goal of this book is to show how to take results of large dynamic finite element models and build small Matlab state space dynamic mechanical models for use in control system models.
The diagram in Figure [1](#figure--fig:hatch00-modal-analysis-flowchart) shows the methodology for analyzing a lightly damped structure using normal modes.
1. deriving the undamped equations of motion in physical coordinates
2. solving the eigenvalue problem, yielding eigenvalues (natural frequencies) and eigenvectors (mode shapes).
These first 2 steps are usually performed using a Finite Element software
3. transform the model from physical coordinate system to modal or principal coordinate system using the eigenvector matrix.
In the modal coordinate system, the original undamped **coupled** equations of motion are transform to the same number of undamped **uncoupled** equations representing the motion of a particular mode of vibration of the system.
4. proportional damping is applied
5. the uncoupled equations are easily solved (as each equation corresponds to a sdof system).
The desired responses are then back-transformed into the physical coordinate system using the eigenvector matrix.
Because finite element models usually have a very large number of states, an important step is the reduction of the number of states while still providing correct responses for the forcing function input and desired output points.
- compute the eigenvalues and eigenvectors (as many as dof in the model)
- the goal is then the reduce the size of the model while still maintaining the desired input/output relationships.
This is done in two steps:
1. reduce the number of dof of the model from the original set to a new set which includes only those dof where forces are applied and where responses are desired
2. reduce the number of modes of vibration used for the solution.
For SISO system, the modes are ranked based on the relative importance to the overall response.
For MIMO system, controllability and observability gramians of the modes are estimated.
Tables [3](#figure--fig:hatch00-n-dof-zeros), [2](#table--tab:notations-eigen-vectors-values) and [3](#table--tab:notations-stiffness-mass) summarize the notations of this document.
The origin and influence of poles are clear: they represent the resonant frequencies of the system, and for each resonance frequency, a mode shape can be defined to describe the motion at that frequency.
We here which to give an intuitive understanding for **when to expect zeros in SISO mechanical systems** and **how to predict the frequencies at which they will occur**.
(<ahref="#citeproc_bib_item_2">Miu 1993</a>) shows that the zeros of any particular transfer function are the poles of the constrained system to the left and/or right of the system defined by constraining the one or two dof's defining the transfer function.
The modal method allows one to replace the n-coupled differential equations with n-uncoupled equations, where each uncoupled equation represents the motion of the system for that mode of vibration.
The overall response of the system is then reconstructed as a superposition of the responses of the different modes of the system.
</div>
Heavily damped structures or structures which explicit damping elements, such as dashpots, result in complex modes and require state space solution techniques using the original coupled equations of motion.
Thus, the present methods only works for lightly damped structures.
Summarizing the modal analysis method of analyzing linear mechanical systems and the benefits derived:
1.**Solve the undamped eigenvalue problem**, which identifies the resonant frequencies and mode shapes (eigenvalues and eigenvectors), useful in themselves for understanding basic motions of the system
2.**Use the eigenvectors to uncoupled or diagonalize the original set of coupled equations**, allowing the solution of n-uncoupled sdof problems instead of solving a set of n-coupled equations
3.**Calculate the contribution of each mode to the overall response**.
This also allows one to reduce the size of the problem by eliminating modes that cannot be excited and/or modes that have not output at the desired dof.
Also, high frequency modes that have little contribution to the system at lower frequencies can be eliminated or approximately accounted for, further reducing the size of the system to be analyzed
4.**Write the system matrix** \\(\bm{A}\\), by inspection.
**Assemble the input and output matrices**, \\(\bm{B}\\) and \\(\bm{C}\\), using appropriate eigenvector terms.
Frequency domain and forced transient response problems can be solved at this point.
If complete eigenvectors are available, initial condition transient problems can also be solved.
For lightly damped systems, proportional damping can be added, while still allowing the equations to be uncoupled.
Having normal modes means that at certain frequencies all points in the system will vibrate at the same frequency and in phase, i.e., **all points in the system will reach their minimum and maximum displacements at the same point in time**.
</div>
Having normal modes can be expressed as:
\begin{equation}
\bm{z}\_i = \bm{z}\_{mi} \sin(\omega\_i t + \phi\_i) = \bm{z}\_{mi} \text{Im}(e^{j\omega\_i t + \phi\_i}) \label{eq:principal\_mode}
\end{equation}
where:
- \\(\bm{z}\_i\\): vector of displacement for all dof's at the i'th frequency
Virtual interpretation of the eigenvectors are shown in Figures [5](#figure--fig:hatch00-tdof-mode-1), [6](#figure--fig:hatch00-tdof-mode-2) and [7](#figure--fig:hatch00-tdof-mode-3).
### Uncoupling the Equations of Motion {#uncoupling-the-equations-of-motion}
At this point, the system is well defined in terms of natural frequencies and modes of vibration.
If any further information such as transient or frequency response is desired, solving for it would be laborious because the system equations are still coupled.
It is thus useful to **transform the n-coupled second order differential equations to n-uncoupled second order differential equations by transforming from the physical coordinate system to a principal coordinate system**.
In linear algebra terms, the transformation from physical to principal coordinates is known as a **change of basis**.
There are many options for change of basis, but we will show that **when eigenvectors are used for the transformation, the principal coordinate system has a physical meaning: each of the uncoupled sdof systems represents the motion of a specific mode of vibration**.
</div>
The n-uncoupled equations in the principal coordinate system can then be solved for the responses in the principal coordinate system using the well known solutions for the single dof systems.
The n-responses in the principal coordinate system can then be **transformed back** to the physical coordinate system to provide the actual response in physical coordinate.
When \\(i \neq j\\), the term \\((\omega\_i^2 - \omega\_j^2)\\) cannot be equal to zero, meaning that:
\begin{equation}
\bm{z}\_{mj}^T \bm{m} \bm{z}\_{mi} = m\_{ij} = 0
\end{equation}
where \\(m\_{ij}\\) is an off-diagonal term in the mass matrix of the principal coordinate system.
The two eigenvectors \\(\bm{z}\_{mi}\\) and \\(\bm{z}\_{mj}\\) are said to be orthogonal with respect to \\(\bm{m}\\).
For \\(i = j\\), \\((\omega\_i^2 - \omega\_j^2) = 0\\), and thus the product \\(\bm{z}\_{mj}^T \bm{m} \bm{z}\_{mi}\\) can be set equal to any arbitrary constant \\(m\_{ii}\\):
\begin{equation}
\bm{z}\_{mi}^T \bm{m} \bm{z}\_{mi} = m\_{ii}
\end{equation}
This is where various **normalization** techniques for eigenvectors come into play.
The stiffness matrix \\(\bm{k}\\), is normalize the same manner.
Because eigenvectors are only known as **ratios** of displacements, not as absolute magnitudes, we can choose how to normalize them.
#### Normalizing with Respect to Unity {#normalizing-with-respect-to-unity}
One method is to normalize with respect to unity, making the **largest** element in each eigenvector equal to unity by dividing each column by its largest value.
Note that the diagonal terms of the stiffness matrix are the squares of the corresponding three eigenvalues.
The normalized stiffness matrix is known as the **spectral matrix**.
Normalizing with respect to mass results in an identify principal mass matrix and squares of the eigenvalues on the diagonal in the principal stiffness matrix, this normalization technique is thus very useful for the following reason.
Since we know the form of the principal matrices when normalizing with respect to mass, no multiplying of modal matrices is actually required: **the homogeneous principal equations of motion can be written by inspection knowing only the eigenvalues**.
</div>
### Transforming Initial Conditions and Forces {#transforming-initial-conditions-and-forces}
Let's determine how to transform initial conditions and forces to the principal coordinate system.
Then, we can solve for transient and forced responses in the principal coordinate system using the uncoupled equations.
We start with the equation of motion in the physical coordinates:
\begin{equation}
\bm{m} \ddot{\bm{z}} + \bm{k} \bm{z} = \bm{F}
\end{equation}
Pre-multiplying by \\(\bm{z}\_n^T\\) and inserting \\(I = \bm{z}\_n \bm{z}\_n^{-1}\\) gives:
- \\(\bm{m}\_p = \bm{z}\_n^T \bm{m} \bm{z}\_n\\) the diagonal principal mass matrix
- \\(\bm{k}\_p = \bm{z}\_n^T \bm{k} \bm{z}\_n\\) the diagonal principal stiffness matrix
- \\(\ddot{\bm{z}}\_p = \bm{z}\_n^{-1} \ddot{\bm{z}}\\) acceleration vector in principal coordinates
- \\(\bm{z}\_p = \bm{z}\_n^{-1} \bm{z}\\) displacement vector in principal coordinates
- \\(\bm{F}\_p = \bm{z}\_n^T \bm{F}\\) force vector in principal coordinates
</div>
The vectors of initial displacements \\(\bm{z}\_{op}\\) and velocities \\(\dot{\bm{z}}\_{op}\\) in the principal coordinate system can be expressed as:
The variables in physical coordinates are the positions and velocities of the masses.
The variables in principal coordinates are the displacements and velocities of each mode of vibration.
</div>
The equations in the principal coordinate system can easily be solved, since the equations are uncoupled.
We now need to **back transform** the results in the principal coordinate system to the physical coordinate system to get the final answer.
We showed previously that the relationship between physical and principal coordinates is:
\begin{equation}
\bm{z}\_n^{-1} \bm{z} = \bm{z}\_p
\end{equation}
where:
- \\(\bm{z}\_n\\) is the matrix containing the eigenvectors (\\(n \times m\\))
- \\(\bm{z}\\) is the vector of physical coordinates (\\(n \times 1\\))
- \\(\bm{z}\_p\\) is the vector of the principal coordinates (\\(m \times 1\\))
And thus, the displacement vector in physical coordinates is obtained by pre-multiplying the vector of displacement in principal coordinates by the normalized modal matrix \\(\bm{z}\_n\\):
\begin{equation}
\bm{z} = \bm{z}\_n \bm{z}\_p
\end{equation}
### Reducing the model size when only selection degrees of freedom are required {#reducing-the-model-size-when-only-selection-degrees-of-freedom-are-required}
The reduction of the number of degrees of freedom is one of the key steps in order to reduce the size of models derived from large finite element simulations.
This can be done as only portions of the eigenvector matrix are needed when only defined dof's have forces applied and other dof's are needed for output.
Let's first examine the force transformation from physical to principal coordinates:
If force is only to be applied on mass 1 (\\(F\_2 = F\_3 = 0\\)), then only the first row of the modal matrix is required to transform the force in physical coordinates to the force in principal coordinates.
Let's now examine the displacement transformation from principal to physical coordinates:
And thus, if we are only interested in the physical displacement of the mass 2 (\\(z\_2 = z\_{n21} z\_{p1} + z\_{n22} z\_{p2} + z\_{n23} z\_{p3}\\)), only the second row of the modal matrix is required to transform the three displacements \\(z\_{p1}\\), \\(z\_{p2}\\), \\(z\_{p3}\\) in principal coordinates to \\(z\_2\\).
**Only the rows of the modal matrix that correspond to degrees of freedom to which forces are applied and/or for which displacements are desired are required to complete the model.**
</div>
### Damping in Systems with Principal Modes {#damping-in-systems-with-principal-modes}
If a mechanical system is designed with a specific viscous damping element, for example a dashpot, then that element can be added to the system as a viscous damper.
The resulting system is linear, but probably does not exhibit normal modes.
This leads to the inability to diagonalize and uncouple the equations of motion.
Damping in typical structures arises from hysteresis losses in the materials as they are strained, in some cases from viscous losses due to structure/fluid interaction but more importantly form relative motion at the interfaces and boundaries where different parts are attached or grounded.
Unless a specific damping element is used in a structural design, most structures have damping which varies from mode to mode and will be in the range of 0.05% to 2% of critical damping.
#### Conditions Necessary for Existence of Principal Modes in Damped System {#conditions-necessary-for-existence-of-principal-modes-in-damped-system}
With a conservative (undamped) system, normal modes of vibration will exist.
In order to have normal modes in a damped system, the modes shapes must be the same as for the undamped case, and the various parts of the system must pass through their minimum and maximum positions at the same instant in time.
A sufficient condition for the existence of damped normal modes is that the **damping matrix be a linear combination of the mass of stiffness matrices**:
Viscous damping in each mode is taken to be an arbitrary percentage \\(\xi\\) of the critical damping \\(c\_{cr}\\):
\begin{equation}
c = \frac{\xi}{2 \sqrt{k m}} = \xi \cdot c\_{cr}
\end{equation}
#### Proportional to Stiffness Matrix: "Relative" Damping {#proportional-to-stiffness-matrix-relative-damping}
Recognizing that the higher modes of vibration damp out quickly, "relative" damping yields damping in proportion to frequencies in normal modes:
\begin{equation}
\xi\_i = \frac{a + b\omega\_i^2}{2 \omega\_i}
\end{equation}
If a value of \\(\xi\_1\\) for the first mode, is assumed:
\begin{equation}
b = \frac{2 \xi\_1}{\omega\_1}
\end{equation}
The value of the damping for any other mode is:
\begin{equation}
\xi\_i = \xi\_1 \frac{\omega\_i}{\omega\_1}
\end{equation}
#### Proportional to Mass Matrix: "Absolute" Damping {#proportional-to-mass-matrix-absolute-damping}
Absolute damping is based on making \\(b = 0\\), in which case the percentage of critical damping is inversely proportional to the natural frequency of each mode.
## Modal Analysis: State Space Form {#modal-analysis-state-space-form}
## Frequency Response: Modal Form {#frequency-response-modal-form}
### Review from Previous Results {#review-from-previous-results}
Since the problem we are solving is frequency response, or finding the steady state motion of each mass as a function of frequency and of applied forces, initial conditions are not required.
From previous analysis, we know the eigenvalues and eigenvectors normalized with respect to mass:
Knowing that in principal coordinates the mass matrix is the identify matrix and the stiffness matrix is a diagonal matrix with the squares of the respect eigenvalues as terms, we can write the matrices by **inspection**:
### Back-Transforming Mode Contributions to Transfer Functions in Physical Coordinates {#back-transforming-mode-contributions-to-transfer-functions-in-physical-coordinates}
The transfer functions in principal coordinates can be back-transformed to physical coordinates.
This allows one to see the contributions of each mode, where \\(z\_{ij}\\) is the physical displacement at dof i due to a force at dof j:
### Forcing Function Combinations to Excite Single Mode {#forcing-function-combinations-to-excite-single-mode}
It is instructive to see what types of forcing function combinations will excite each of the three modes separately.
From the definition of normal modes, we known that if the system is started from initial displacement conditions that match one of the normal modes, the system will respond at that only mode.
An analogous situation exists for combinations of forcing functions.
The forces transform in the principal coordinates using:
Thus, if \\(\bm{F}\\) is aligned with \\(\bm{z}\_{ni}\\) (the i'th normalized eigenvector), then \\(\bm{F}\_p\\) will be null except for its i'th term and only the i'th mode will be excited.
</div>
### How Modes Combine to Create Transfer Functions {#how-modes-combine-to-create-transfer-functions}
Any transfer function derived from the modal analysis is an additive combination of sdof systems.
Each single degree of freedom system has a gain determined by the appropriate eigenvector entries and a resonant frequency given by the appropriate eigenvalue.
It can be shown that for a general system with \\(m\\) undamped modes:
Equations <eq:general_add_tf> and <eq:general_add_tf_damp> shows that in general every transfer function is made up of **additive combinations of single degree of freedom systems**, with each system having its DC gain determined by the appropriate eigenvector entry product divided by the square of the eigenvalue, \\(z\_{nji} z\_{nki}/\omega\_i^2\\), and with resonant frequency defined by the eigenvalue \\(\omega\_i\\).
{{<figuresrc="/ox-hugo/hatch00_z11_tf.png"caption="<span class=\"figure-number\">Figure 10: </span>Mode contributions to the transfer function from \\(F\_1\\) to \\(z\_1\\)">}}
The zeros for SISO transfer functions are the roots of the numerator, however, from modal analysis we can see that the zeros arise when modes combine with appropriate phase such that the resulting motion is null.
## SISO State Space Matlab Model from ANSYS Model {#siso-state-space-matlab-model-from-ansys-model}
In this section is developed a SISO state space Matlab model from an ANSYS cantilever beam model as shown in Figure [11](#figure--fig:hatch00-cantilever-beam).
{{<figuresrc="/ox-hugo/hatch00_cantilever_beam.png"caption="<span class=\"figure-number\">Figure 11: </span>Cantilever beam with forcing function at midpoint">}}
ANSYS has a number of different eigenvalue extraction techniques, but for most problems only two methods are commonly used:
- The first method, **Block Lanczos**, is the fastest and calculates all the eigenvalues or eigenvectors in a specific frequency range.
This method is usually preferred as most practical models require knowledge of the modes from DC through a specified higher frequency.
- The second method, Reduced, performs a **Guyan reduction** on the model to reduce its size, then calculates all the eigenvalues for the reduced model.
Obtaining eigenvector components for the reduced degrees of freedom requires an additional calculation step in ANSYS.
This method is still useful to obtain a "super-element".
We can choose to use ANSYS to output only the eigenvectors for nodes of interest or we can output the complete modal matrix and choose the appropriate rows of data within Matlab.
### Matlab State Space Model from ANSYS Eigenvalue Run {#matlab-state-space-model-from-ansys-eigenvalue-run}
In order to obtain a State Space model of reasonable order, we need to rank the relative importance of the contributions of each of the individual modes.
To do so, we can use a **ranking of DC gains**.
Once the modes are ranked, the most important can be selected for use, and modes with small DC gains are eliminated from the model.
The DC gain contributions of the eliminated modes are not included in the overall DC gain, so there is error in the low frequency gain.
In order to eliminate this error, the Matlab function `modred` is introduced.
Using `modred` is analogous to using Guyan reduction to reduce some less important degrees of freedom.
We will discuss in this section two methods of sorting, one which is applicable for models with the same value of damping for all modes \\(\xi\_i = \xi\\) ("uniform" damping), and another which is applicable for models with different damping values for each mode ("non-uniform" damping).
The general equation for the overall transfer function of undamped and damped systems are:
And we see that the peak gain for a mode is the DC gain of the same mode divided by \\(2 \xi\_i\\).
If the same value of \\(\xi\\) is used for all modes, then all the DC gain terms are divided by the same \\(2 \xi\\) terms and the relative amplitude of the DC gains and peak gains are the same, so there is no difference between sorting a uniform damping model using DC gain or peak gain.
However, if the modes have different damping the relationship between the DC gain and peak gain for all the modes is not a constant value and peak gain must be used to rank modes for importance.
The Matlab command `modred` can be used for reducing models while retaining the overall system DC gain.
The `matchdc` gain option for the function `modred` reduces defined states by setting the derivatives of the states to be eliminated to zero, then solving for the remaining states.
The method essentially sets up the eliminated states to be "infinitely fast" and is analogous to Guyan reduction in that the low frequency effects of the eliminated states are included in the remaining states.
In thus case, `modred` produces a new state space model with a non-null \\(D\\) matrix meaning that the reduced system will have a flat response at high frequency which may not be desirable.
The `truncate` method does not account for the DC gains of the unused modes, which can result in error in the low frequency portion of the frequency response.
However, the `truncate` method has the advantage that it does not exhibit the unusual high frequency direct transmission matrix related behavior of the `matchdc` method.
If sorting of DC gain values is performed prior to the `truncate` operation, the system DC gain error may be acceptable while maintaining better high frequency performance.
## Ground Acceleration Matlab Model From ANSYS Model {#ground-acceleration-matlab-model-from-ansys-model}
In this section we wish to extract a SISO state space model from a Finite Element model representing a Disk Drive Actuator (Figure [12](#figure--fig:hatch00-disk-drive-siso-model)).
The primary motion of the actuator is rotation about the pivot bearing, therefore the final model has the coordinate system transformed from a Cartesian x,y,z coordinate system to a Cylindrical \\(r\\), \\(\theta\\) and \\(z\\) system, with the two origins coincident (Figure [13](#figure--fig:hatch00-disk-drive-nodes-reduced-model)).
{{<figuresrc="/ox-hugo/hatch00_disk_drive_nodes_reduced_model.png"caption="<span class=\"figure-number\">Figure 13: </span>Nodes used for reduced Matlab model. Shown with partial finite element mesh at coil">}}
The four nodes 24061, 24066, 24082 and 24087 are located in the center of the coil in the z direction and are used for simulating the VCM force.
The arrows at the nodes indicate the direction of forces.
Nodes 22 and 10022 are the nodes for the top and bottom heads and are used for measurement.
### Ansys Actuator/Suspension Model Results {#ansys-actuator-suspension-model-results}
A recommended sequence for analyzing dynamic finite element models is:
1. Plot resonant frequencies versus mode numbers to get a feel for the frequency range.
See if there are any significant jumps in frequency between modes which can indicate the system transitioning from one type of characteristic motion to another.
2. Plot frequency responses to define which modes couple into the response.
3. Plot and animate the mode shapes that contribute to the response, identifying modes that couple into motions in directions of interest and those that do not.
Visually get a sense of how the geometry of the structure affects the modes.
4. Run parameter studies to understand the sensitivity of critical modes to design variables: dimensions, tolerances, material properties, etc.
### Ansys Output Example Listing {#ansys-output-example-listing}
A small section of the exported `.eig` file from ANSYS is shown bellow..
- The eigenvector for the mode in all dofs for all nodes
### Inputs and Outputs definition {#inputs-and-outputs-definition}
A unity force is applied at the coil, and evenly distributed among the four nodes.
This model is still considered as a **Single Input** model because the same force is applied to all four coil nodes, requiring only a single column vector for the input matrix \\(\bm{B}\\).
To compute the DC gain and peak gain in that case, we define a **composite** forcing function which consists of the force applied to each node times the eigenvector value for that node, \\(\bm{F}\_p^T \bm{x}\_n\\).
The dimensions of this operation are \\((1 \times n) \times (n \times m) = (1 \times m)\\), so we have a composite force vector for each mode.
### Building Full State Space Matrices {#building-full-state-space-matrices}
From Ansys, we have the eigenvalues \\(\omega\_i\\) and eigenvectors \\(\bm{z}\\).
In this chapter another method of reducing models, “balanced reduction”, will be introduced and compared with the DC and peak gain ranking methods.
This method uses the concepts of controllability and observability, commonly referenced in the control community.
One issue with balanced reduction is that we lose the ability to directly identify individual modes in the reduced system model.
After balanced reduction one needs to examine the system matrix to identify which modes are included, while the dc and peak gain ranking techniques retain the identities of the individual modes.
Unlike SISO models, which can be easily ranked using simple DC and peak gain techniques, **MIMO models will require the balanced reduction method because it easily handles the problem of ranking multiple inputs and outputs**.
### Reviewing DC gain raking {#reviewing-dc-gain-raking}
So far we have used DC or peak gains of the individual modes to rank the importance of including each mode in the reduced system.
For any mode, if the degree of freedom associated with the applied force has a zero value, then the force applied at the degree of freedom cannot excite that mode, so the DC and peak gains will also be zero and the mode be eliminated.
Similarly, if the degree of freedom associated with the output has a zero value, then no matter how much force is applied to that mode, there will be no output, and the mode can be eliminated.
A mode which cannot be excited by the applied force is said to be **uncontrollable**, and a mode which has no output in the desired direction is said to be **unobservable**.
the following definitions of controllability hold:
1. If there is an input \\(u\\) that can move the system from some arbitrary state \\(\bm{x}\_1\\) to another arbitrary state \\(\bm{x}\_2\\) in a finite time, then the system is controllable
2. A controllability matrix \\(\bm{\mathcal{C}}\\) can be formed as:
\bm{W}\_{c} = \int\_{0}^{\infty} e^{\bm{A}\tau} \bm{B} \bm{B}^{T} e^{{\bm{A}^{T} \tau} d \tau}
\end{equation}
If the solution \\(\bm{W}\_c\\) is non-singular, then the system is controllable.
Diagonal elements of the controllability gramian give information about the **relative controllability of the different modes** and can be used in a manner similar to our use of DC gains to rank the relative controllability of individual modes.
Gramians exists only for system that have all their poles with negative real parts.
Thus, if a system has rigid body modes, it should first be partitioned into rigid body dynamics and oscillatory dynamics and then the controllability gramian can be computed on the oscillatory partition.
A similar set of definitions can be made for observability:
1. If the initial state \\(\bm{x}\_0\\) of a system can be inferred from knowledge of the input \\(u\\) and output \\(\bm{y}\\) over a finite period of time \\((0, t)\\), then the system is said to be observable.
2. A observability matrix \\(\bm{\mathcal{O}}\\) can be formed as:
If the solution \\(\bm{W}\_o\\) is non-singular, then the system is observable.
Diagonal elements of the observability gramian give information about the **relative observability of the different modes** and can be used in a manner similar to our use of DC gains to rank the relative observability of individual modes.
### Ranking Using Controllability/Observability {#ranking-using-controllability-observability}
We could use the controllability curve to rank the states for controllability and eliminate those states with low controllability.
Alternately, we could use the observability curve to rank the states for observability and then eliminate states with low observability.
The problem with this approach is that the joint controllability/observability is not taken in account.
To take into account both the observability and the controllability, we use the **Balance Reduction**.
### Balanced Reduction {#balanced-reduction}
The balanced reduction is used to create a system with **identical diagonal controllability and observability Gramians**.
Since the two gramians are equal, either the diagonal of the controllability gramian or the diagonal of the observability gramian can be used to rank states for elimination.
The diagonal terms of the joint gramian are squares of the Hankel singular values of the system.
The Hankel matrix is the product of the controllability and observability gramians.
Because the controllability and observability gramians are identical, there is no ambiguity in deciding whether the most controllable or the most observable states should be chosen.
The **states to be kept are the states with the largest diagonal terms**.
## MIMO Two Stage Actuator Model {#mimo-two-stage-actuator-model}
Instead of applying voltage as the input into the piezo elements, we will assume that we have calculated an equivalent set of forces which can be applied at the ends of the element that will replicate the voltage force function.
In this model, we will be applying forces to multiple nodes at the ends of both piezo elements.
</div>
Since the same forces are being applied to both piezo elements, they represent the second input to the MIMO system, the first input being the coil force.
### Ansys Model Description {#ansys-model-description}
{{<figuresrc="/ox-hugo/hatch00_disk_drive_mimo_ansys.png"caption="<span class=\"figure-number\">Figure 15: </span>Nodes used for reduced Matlab model, shown with partial mesh at coil and piezo element">}}
For a SISO system, we can rank the relative importance of modes using two methods, by using DC or peak gains and by using balancing.
However, **for a MIMO system, balancing is the only practical solution**.
If we were using DC gains to rank the modes, we would have to compute DC gains for the four combinations possible of the two inputs and two inputs.
We would then see that the DC gains of the modes are not ranked to same way for different sets of input/output.
Especially, the modes important for the VCM are not the same than for the piezo.
If one were to choose a single ranking for the model which would take into account both inputs and both outputs, it is difficult to see how to do it given the DC gain rankings.
Thus the necessity of balanced reduction for MIMO models.
#### Balancing Reduction {#balancing-reduction}
Balancing the system involves calculating gramians, which are only defined for negative definite systems.
This requires separating the rigid body mode from the oscillatory modes and balancing the oscillatory modes.
The system matrices are partitioned and a model of only oscillatory modes is created and balanced.
Plotting the diagonal gramian terms reveals the relative important of the states.
The complete system is rebuilt by augmenting the rigid body mode with the reduced oscillatory modes.
{{<figuresrc="/ox-hugo/hatch00_z13_tf.png"caption="<span class=\"figure-number\">Figure 16: </span>Mode contributions to the transfer function from \\(F\_1\\) to \\(z\_3\\)">}}
{{<figuresrc="/ox-hugo/hatch00_z11_tf.png"caption="<span class=\"figure-number\">Figure 17: </span>Mode contributions to the transfer function from \\(F\_1\\) to \\(z\_1\\)">}}
{{<figuresrc="/ox-hugo/hatch00_cant_beam_gramian_balanced.png"caption="<span class=\"figure-number\">Figure 22: </span>Sorted values of the Gramian of the balanced realization">}}
