%Date: Tue, 18 Jul 2000 14:53:41 -0500 (CDT)
%Return-Path: <cha@adder.eng.hmc.edu>
%From: Philip Cha <cha@eng.hmc.edu>
%Subject: no subject (file transmission)
%To: depillis@mcs.anl.gov
%Date: Mon, 17 Jul 100 21:29:21 -0700 (PDT)
%  modified from add.tex

%  combined addnew.tex and add2.tex

%7/17/00

%\bf{h}}ginal 3/18/99, P.D. Cha
% Updated  4/7/99, L.G. de Pillis
% Updated  4/16/99, L.G. de Pillis
% Updated  5/3/99, P. D. Cha
% Updated 5/19/99, L.G. de Pillis
% Updated 5/25/99, L.G. de Pillis
% Updated 5/28/99, P. D. Cha
% Updated 6/21/99, P. D. Cha
% Updated 6/23/99, L.G. de Pillis

% Updated 7/28/99, P. D. Cha
% Updated 8/5/99, P. D. Cha
% Updated 8/5/99, L.G. de Pillis
% Updated 8/5/99, P. D. Cha




%\relax 
%\documentstyle[11pt,doublespace]{article}
%\documentstyle[11pt]{article}
%LdeP Modifications to fit IJNME LaTeX reqs
\documentclass{nmeauth}
%\usepackage{epsfig}
%\usepackage{psfig}

\RequirePackage[dvips]{graphics}


% Would like to include AMSMATH

%\columnsep 0.375in %LdeP
%\topmargin 0.5625in %LdeP
%\linespace 0.2in
%\headheight 0.0in %LdeP
%\headsep 0.0in %LdeP
%\textheight 8.00in %LdeP
%\textwidth 6.5in %LdeP
%\footskip 0.30in %LdeP
%\footheight 0.30in %LdeP
%\renewcommand{\baselinestretch}{2} %LdeP
%\parskip 0.15in %LdeP
%\oddsidemargin 0.25in %LdeP
%\evensidemargin 0.25in %LdeP
%\marginparsep 0.0in %LdeP
%\marginparwidth 0.0in %LdeP

\def\b{\begin{equation}}
\def\e{\end{equation}}
\def\x{\underline x}
\def\X{\underline X}
\def\q{\underline q}
\def\al{\underline\alpha}
\def\y{\underline y}
\def\Y{\underline Y}
\def\h{\hspace{0.3in}}
%\def\y{\bf y}
\def\l{\lambda}
\def\g{\underline\gamma}
\def\D{\underline D}
\def\ep{\epsilon}
\def\v{\vspace{0.25in}}
%\def\baselinestretch{1.0}
\def\baselinestretch{2}

\def\cond{{\rm cond}}

%\pagestyle{empty}

\begin{document}

\NME{1}{45}{0}{0}{00}

\runningheads{P.D. Cha and L.G. de Pillis}{Model Updating with Known Masses}

\received{August 1999}
\revised{August 2000}
\accepted{August 2000}

%\voffset=0.5in

\title{Model Updating by Adding Known Masses}

\author{P.D. Cha\affil{1} and L.G. de Pillis\affil{2}\comma\corrauth}

\address{\affilnum{1}\ Department of Engineering,
Harvey Mudd College, Claremont, CA 91711 \\
\affilnum{2}\ Argonne National Laboratory, Division of Mathematics
and Computer Science, Argonne, IL 60439 and
Department of Mathematics, Harvey Mudd College, Claremont, CA
91711}

\corraddr{Department of Mathematics, Harvey Mudd College,
Claremont, CA 91711, depillis@hmc.edu}


%\centerline{\bf MODEL UPDATING BY ADDING KNOWN MASSES}% AND STIFFNESSES,}
%\centerline{\bf PART I:  THEORETICAL FORMULATION}
%\vspace{1.5in}
%\centerline{Philip D. Cha\footnote{Associate Professor of
%Engineering, Harvey Mudd College, Claremont, CA 91711, USA}}
%\centerline{Lisette G.~de Pillis\footnote{Associate Professor of
%Mathematics, Harvey Mudd College, Claremont, CA 91711, USA}}




%\vspace{0.5in}
%\centerline{Department of Engineering}
%\centerline{Harvey Mudd College}
%\centerline{Claremont, California 91711}


%\vspace{1.25in}
%\noindent\newline{Number of pages:  52}
%%\newline{Number of tables: 5}
%%\newline{Number of figures: 1}
%\newline{Number of appendix: 1}
%\newline{Contact person:  Philip D. Cha (see postal address above)}
%

%\newpage 

\begin{abstract}


New approaches are developed that use 
measured data to adjust the analytical mass
and stiffness matrices of a system so that the
agreement between the analytical modes of vibration
and the modal survey is improved.
%such that the corrected model will have modes
%of vibration that are closer to the modal survey.
By adding known masses to
the structure of interest, measuring the
modes of vibration of this mass-modified system, and finally
using this set of new data in conjunction with the initial modal
survey, the analytical mass matrix of the structure
can be corrected, after which the analytical stiffness matrix
can be readily updated.  Manipulating the correction
matrices into vector forms,
the connectivity information can be enforced,
thereby preserving the physical configuration of the
system and reducing the sizes of the least squares 
problems that need to be solved.  Solution techniques
for updating the system matrices are introduced, and the
numerical issues associated with solving overdetermined and
underdetermined least squares problems are investigated. 
The effects of round-off errors are also studied, and
heuristic criteria are given for determining the minimum number of modes that
need to be measured in order to ensure
sufficiently accurate updated mass and stiffness matrices.
Numerical experiments are presented to validate the proposed
model updating techniques, to illustrate
the effects of the number of measured modes on the
quality of the updated model, to show how the magnitudes and
locations of the added masses influence the
updated matrices, 
and to highlight the numerical issues discussed in this paper.

\end{abstract}

\keywords{Model updating, Added masses, Least squares}


\section{Introduction}

Highly accurate and detailed analytical models 
are required
to analyze and predict the dynamical behavior of 
complex structures during analysis and design.  
With the proliferation of digital 
computers, new methods of analysis have 
been developed, in particular through the method of finite elements.  
Once the finite element model of a physical system
is constructed, its accuracy is often tested
by comparing its analytical modes of vibration 
(or natural frequencies and mode
shapes) with those obtained from the physical system
during a modal survey.  If the agreement
between the two is good, then more 
credence is given to the analytical model, and it
can be used with confidence for future analysis. 
If the correlation between the two is poor,
then assuming that the experimental measurements are 
correct, the analytical model must be 
adjusted so that the agreement between predictions 
and test results is improved.
The updated model may then be considered 
to be a better dynamical representation of the 
structure than the initial analytical model. 
The updated model can subsequently be used
with reasonable accuracy to assess the stability and 
control characteristics and to
predict the dynamical responses of the structure.
The above process of correcting the system
matrices is known as {\it model updating}.


In recent years various methods have been developed
to improve the quality
of the analytical finite element models using test data.  
A detailed discussion of every approach
is beyond the scope of this paper, and interested 
readers are referred to the recent survey paper by
Mottershead and Friswell \cite{Mottershead}.  
In the following paragraphs, some
commonly used model updating techniques will
be briefly reviewed.

%
%
Berman \cite{Berman} proposed an updating scheme
based on the Lagrange multipliers formalism
that uses the measured
mode shapes to correct the mass 
matrix of a structure.  This updating algorithm
identifies a set of minimum changes
in the analytical mass matrix so that the 
measured modes are orthogonal to the updated
mass matrix of the system.
Using essentially the method first introduced by
Baruch and Bar Itzhack \cite{Baruch}, Wei \cite{Wei} 
developed an optimal method to
update the stiffness matrix of a structure.  He
also employed the Lagrange multipliers approach to
adjust the stiffness matrix, subjected to 
the constraints of satisfying the
generalized eigenvalue problem, the orthogonality 
condition of the measured mode shapes
and the symmetry property of the stiffness
matrix.
%

%
The Lagrange multipliers approaches to update the system 
matrices return fully populated (or {\em dense}) mass and stiffness
matrices that may not bear any resemblance to the physical
system being analyzed.  To preserve the physical load
paths of the original analytical model,
Kabe \cite{Kabe} assumed the analytical mass matrix to be correct
and incorporated the readily available structural connectivity
information in addition to the test data
to optimally adjust the stiffness matrix.  
The adjustments he performed ensure that 
zero and nonzero elements 
of the analytical model are preserved, 
and the adjusted model exactly reproduces the modes used 
in the identification.  He also utilized 
a Lagrange multipliers technique, so that 
the percentage change to each
stiffness element is minimized.
While Kabe's approach to updating the
stiffness matrix is straightforward, 
the assumption that the actual
mass matrix is identical to the analytical mass matrix 
remains debatable \cite{Berman2}.
%, since in practice
%there is significantly greater success using
%finite elements during static analyses (which use
%only the stiffness matrix) than during
%dynamic analyses (which use both the mass and stiffness matrix).  
%

Using an approach based on matrix
perturbation theory, Chen, Kuo and Garba \cite{ChenKuo} found the
correction mass and stiffness matrices
by enforcing the orthogonality conditions of the
measured mode shapes with respect to the system
matrices.   Like the schemes proposed by Berman and Wei, however, 
the updating algorithms also return fully
populated  mass and stiffness
matrices, thus failing to preserve the physical
connectivity of the system.  In addition, because the approach
outlined in \cite{ChenKuo} is based on perturbation theory, the updating
algorithm can only be applied when the deviations
of the actual parameters from the analytical values are small.
%Otherwise, the underlying assumptions allowing for asymptotic expansions
%to be used in perturbation theory will be violated, and the
%perturbation updating scheme will yield erroneous results.
Finally, the derivation carried out by Chen {\it et al.}
requires that the measured modal matrix, $[X]$, be properly
normalized with respect with the actual mass and stiffness
matrices, $[M]$ and $[K]$, of the system, such that
it satisfies $[X]^T[M][X]=[I]$ and $[X]^T[K][X]=[\Lambda]$.  Because the objective of
model updating is to correct the system matrices, $[M]$ and $[K]$
are not known {\it a priori}.  Thus, the orthogonality constraints
cannot be enforced, and their updating algorithm cannot be applied in
practice.
%
%
%

In this paper, new model updating schemes are introduced that adjust
the system mass and stiffness matrices from an 
incomplete set of measured modes.
The underlying principle of the proposed scheme is to
add known masses to the physical structure,
measure the modes of vibration of the new system, and then
use the new set of measurements in conjunction with
the original set of experimental data to
correct the mass and stiffness matrices of the actual
structure.  

%Nalitolela, Penny and Friswell [8] 
%added masses and stiffnesses to the system
%in an attempt to update the structural parameters.  However,
%they did so to overcome the problems of incomplete modes 
%and the mismatch between the degrees-of-freedom of the analytical
%model and the number of measured coordinates during testing.
%Here, the modes of vibration of the mass-modified
%system are used in conjunction with the modal
%characteristics of the initial system to directly update the
%analytical system matrices. 



\section{Proposed Model Updating Algorithm}

Consider the analytical model of a given structure,
with $N$ degrees-of-freedom, whose modes of vibration
are given by the solutions of the following generalized
eigenvalue problem: 
\begin{equation} \label{eigan}
[K_o][X_o]=[M_o][X_o][\Lambda_o]
\end{equation}
where $[M_o]$ and $[K_o]$ are the symmetric analytical
mass and stiffness matrices of the system,
$[X_o]$ is the $N\times N$ modal matrix (whose
columns correspond to the eigenvectors or mode shapes) 
of the analytical model, and $[\Lambda_o]$ is an $N\times N$
diagonal matrix whose elements correspond to the 
eigenvalues (or the square of the natural frequencies) of the
analytical model.  A mode of vibration constitutes a given
natural frequency and its corresponding mode shape (or the
square root of an eigenvalue and its corresponding eigenvector).

Experimentally, it is often difficult if not impossible to measure
the same number of modes as the number
of degrees-of-freedom of the analytical model.  
Thus, the measured data are said to be incomplete.
A problem unrelated to that previously described, but
also commonly referred to as ``incomplete'' occurs
when the measured eigenvector contains fewer coordinates
than are available from the analytical model.  In common
with other model updating techniques, 
%this case, 
the measured eigenvectors must first be expanded
before the proposed algorithms can be applied.  Various
mode expansion techniques can be found in references
\cite{Gysin} and \cite{Imregun}.
In this paper, we will assume that all the coordinates of the
eigenvectors can be measured. 
This allows us to focus our attention on the quality of the proposed
updating algorithm, and not confound the resulting updates with
errors introduced by mode shape expansion.  Therefore,
we will reserve the
word ``incomplete'' to mean when the test measurements
contain fewer modes than those of the analytical model.







Regardless of whether the measured modes are complete or incomplete, 
the modes of vibration of the actual system must satisfy
the following generalized eigenvalue problem:
\begin{equation} \label{eig}
[K][X]=[M][X][\Lambda]
\end{equation}
where $[M]$ and $[K]$ are the actual $N\times N$
symmetric mass and
stiffness matrices of the physical system, $[X]$ is 
the measured $N\times N_e$ modal matrix ($N_e$ denotes
the number of measured modes; $N_e \le N$), 
and $[\Lambda]$ is an $N_e\times N_e$ diagonal matrix
whose elements correspond to the measured eigenvalues of the
system.  Knowing $[M_o]$, $[K_o]$, $[X]$ and $[\Lambda]$, the 
objective of model updating is to correct $[M_o]$ and $[K_o]$
so that the new analytical system matrices yield
modes of vibration that are closer to the measured data than they
were initially.
%The approaches we implement to update the mass and stiffness matrices
%are outlined below.


\subsection{Mass Updating}\label{massupdating}

To update the mass matrix of the analytical model, 
we add a known mass matrix, 
$[M_a]$, to the physical structure, at locations coincident
with the nodes of the finite element model in
order to preserve the size of the initial analytical
model.  The resulting
system satisfies
\begin{equation} \label{18}
[K][X_a]=\left( [M]+[M_{a}] \right)[X_a][\Lambda_a]
\end{equation}
where $[X_a]$ corresponds to the $N\times N_e$ measured modal
matrix of the new system, and $[\Lambda_a]$ is an $N_e\times N_e$
diagonal matrix, whose elements are the measured eigenvalues 
of this mass-modified system.  
It is assumed that the added masses do not alter the stiffnesses
of the system significantly.
Taking the transpose of Eq.~(\ref{eig}) and postmultiplying the
resulting matrix equation by $[X_a]$, we get
\begin{equation} \label{aa}
[X]^T[K][X_a]=[\Lambda][X]^T[M][X_a]
\end{equation}
Premultiplying Eq.~(\ref{18}) by $[X]^T$, we have
\begin{equation} \label{bb}
[X]^T[K][X_a]=[X]^T\left( [M]+[M_{a}] \right)[X_a][\Lambda_a]
\end{equation}
Subtracting Eq.~(\ref{bb}) from (\ref{aa}), 
%Equating the right hand sides of Eqs.~(\ref{aa}) and (\ref{bb}), 
we obtain
%\begin{equation} 
%[\Lambda][X]^T[M][X_a]=[X]^T\left( [M]+[M_{a}] \right)[X_a][\Lambda_a]
%\e
%or
\begin{equation} \label{c}
[\Lambda][X]^T[M][X_a]-[X]^T[M][X_a][\Lambda_a]=[Q]
\end{equation}
where
\begin{equation} 
[Q]=[X]^T [M_{a}] [X_a][\Lambda_a]
\end{equation}
Let
\begin{equation} \label{23}
[P]=[X]^T[M][X_a]
\end{equation}
then Eq.~(\ref{c}) simplifies to
\begin{equation} \label{dd}
[\Lambda][P]-[P][\Lambda_a]=[Q]
\end{equation}
Because both $[\Lambda]$ and $[\Lambda_a]$
are diagonal matrices, Eq.~(\ref{dd}) can be easily
expanded so that its $(i,j)$th element
yields
\begin{equation} \label{26}
(\lambda_i - \lambda_{aj} )P_{ij}=Q_{ij}
\end{equation}
where $\lambda_{aj}$ is the $j$th measured eigenvalue
of the mass-modified system and 
$i,j=1,\ldots,N_e$.  Assuming the
$N_e$ measured eigenvalues of the original and the
mass-modified systems do not coincide,
we can solve for all of the unknowns $P_{ij}$
and then construct matrix $[P]$.  If any two 
measured eigenvalues of the 
original and the mass-modified systems coincide,
we simply change the added masses or their locations
to make the eigenvalues distinct.  
%See section 5.4
See Section \ref{effects}
for further discussion on the selection
of the number, placement and magnitude of the added masses.
%We discuss the selection of the number, placement and magnitude
%of the added masses in the section on numerical experiments.







%LdeP QUESTION
%******** QUESTION: HOW CAN WE GUARANTEE THAT $(\lambda_i - \lambda_{aj} )<>0$
%for all pairs $(i,j)$?

The actual and the analytical mass matrices of the system
are related as follows
\begin{equation} 
[M] = [M_o] + [\delta M]
\end{equation}
where $[\delta M]$ represents the correction to the analytical 
mass matrix.  Rewriting Eq.~(\ref{23}), we obtain
\begin{equation} \label{massnew1}
[X]^T[\delta M][X_a]=[P]-[X]^T[M_o][X_a]
\end{equation}
Because $[X]$ and $[X_a]$ are both rectangular
matrices (assuming $N_e<N$), they have
no inverses.  However, 
Eq.~(\ref{massnew1}) can be rewritten so that $[\delta M]$ appears
as an unknown column vector $\delta\bf{m}$ as follows:
\begin{equation} \label{massnew}
[\,A\,] \, \delta\bf{m}=\bf{r}
\end{equation}
where 
\begin{equation} \label{dm}
\delta{\bf m}=\left[
\delta m_{11} \cdots \delta m_{1N} \;|\;
\delta m_{21} \cdots \delta m_{2N} \;|\; \cdots \;|\;
\delta m_{N1} \cdots
\delta m_{NN} \right]^T
\end{equation}
and 
\begin{equation} 
{\bf r}=\left[
r_{11} \cdots r_{1N_e} \;|\;
r_{21} \cdots r_{2N_e} \;|\; \cdots \;|\;
r_{N_e1} \cdots
r_{N_eN_e} \right]^T
\end{equation}
In Eq.~(\ref{dm}), $\delta m_{ij}$ corresponds
to the $(i,j)$th element of $[\delta M]$.
Matrix $[A]$ is of size $N_e^2\times N^2$,
whose elements can be determined by expanding
the left hand side of Eq.~(\ref{massnew1}); vector $\bf r$ is of length
$N_e^2$, whose components can be obtained by expanding the
right hand side of Eq.~(\ref{massnew1}).
The technique to solve Eq.~(\ref{massnew}) will be discussed 
in Section \ref{solutiontechnique}.
%in a separate section.
We now turn our attention to updating the stiffness matrix
of the system.


\subsection{Stiffness Updating}

From Eqs.~(\ref{aa}) and (\ref{bb}), we get
\begin{equation} \label{aas}
[\Lambda]^{-1}[X]^T[K][X_a]=[X]^T[M][X_a]
\end{equation}
and
\begin{equation} \label{bbs}
[X]^T[K][X_a][\Lambda_a]^{-1}=[X]^T\left( [M]+[M_{a}] \right)[X_a]
\end{equation}
Subtracting Eq.~(\ref{aas}) from Eq.~(\ref{bbs}), we have
\begin{equation} \label{jj}
[X]^T[K][X_a][\Lambda_a]^{-1}-[\Lambda]^{-1}[X]^T[K][X_a]=[S]
\end{equation}
where
\begin{equation} 
[S]=[X]^T[M_a][X_a]
\end{equation}
Defining
\begin{equation} \label{243}
[U]=[X]^T[K][X_a]
\end{equation}
then Eq.~(\ref{jj}) simplifies to
\begin{equation} \label{w3}
[U][\Lambda_a]^{-1}-[\Lambda]^{-1}[U]=[S]
\end{equation}
Because both $[\Lambda]$ and $[\Lambda_a]$ are diagonal,
Eq.~(\ref{w3}) can be easily expanded so that its
$(i,j)$th element yields
\begin{equation} 
\left(\frac{1}{\lambda_{aj}} - \frac{1}{\lambda_i} \right)U_{ij}=S_{ij}
\end{equation}
%
%
%==================================================================
%
%To update the stiffness matrix of the analytical model,
%we add a known stiffness matrix, $[K_b]$, to the actual 
%system, at locations coincident with the nodes of the finite
%element model in order to maintain the size
%of the initial analytical model, such that the 
%stiffness-modified system satisfies
%\begin{equation} \label{28}
%([K]+[K_b])[X_b]=[M][X_b][\Lambda_b]
%\e
%where $[X_b]$ corresponds to the $N\times N_e$ modal
%matrix of the new system, and $[\Lambda_b]$ consists
%of an $N_e\times N_e$  diagonal matrix
%whose elements are the eigenvalues 
%of this stiffness-modified system.
%Both matrices are experimentally determined.
%Here, it is implicitly assumed that the additional
%stiffness matrix does not significantly affect
%the mass matrix of the system.
%Taking the transpose of Eq.~(\ref{eig}) and postmultiplying the
%resulting matrix equation by $[X_b][\Lambda_b]$, we get
%\begin{equation} \label{2aa}
%[X]^T[K][X_b][\Lambda_b]=[\Lambda][X]^T[M][X_b][\Lambda_b]
%\end{equation}
%Premultiplying Eq.~(\ref{28}) by $[\Lambda][X]^T$, we have
%\begin{equation} \label{2bb}
%[\Lambda][X]^T([K]+[K_b])[X_b]=
%[\Lambda][X]^T [M][X_b][\Lambda_b]
%\end{equation}
%Subtracting Eqs.~(\ref{2aa}) and (\ref{2bb}), we obtain
%\begin{equation} \label{2c}
%[X]^T[K][X_b][\Lambda_b]-[\Lambda][X]^T[K][X_b]
%=[S]
%\end{equation}
%where
%\begin{equation} 
%[S]=[\Lambda][X]^T[K_b] [X_b]
%\end{equation}
%Defining 
%\begin{equation} \label{243}
%[U]=[X]^T[K][X_b]
%\end{equation}
%then Eq.~(\ref{2c}) simplifies to
%\begin{equation} \label{2dd}
%[U][\Lambda_b]-[\Lambda][U]=[S]
%\end{equation}
%Because both $[\Lambda]$ and $[\Lambda_b]$
%are diagonal, Eq.~(\ref{2dd}) can be easily
%expanded so that its $(i,j)$th element
%yields
%\begin{equation} \label{29}
%(\lambda_{bj}- \lambda_{i} )U_{ij}=S_{ij}
%\end{equation}
where $i,j=1,\ldots,N_e$.  Assuming that the $N_e$ 
measured eigenvalues of the original and the mass-modified
structures are distinct,   
%LdeP QUESTION
%QUESTION: AGAIN, HOW CAN WE GUARANTEE THAT $(\lambda_{bj}- \lambda_{i} )$
%is NEVER ZERO?
the unknowns $U_{ij}$ can be computed and
matrix $[U]$ can be assembled.
%If any two measured eigenvalues $\lambda_{bj}$ and $\lambda_{i}$
%coincide, we simply change the added stiffnesses or their locations
%to make the eigenvalues distinct.
%The selection of the placement, magnitude and number of added stiffnesses
%can be chosen in a manner similar to that of the added masses.  See the
%companion paper, Part II.

As in the case of the mass matrix, 
the actual and the analytical stiffness matrices
are related by
\begin{equation} 
[K]=[K_o]+[\delta K]
\end{equation}
where $[\delta K]$ represents the matrix of stiffness offsets.
In order to find $[\delta K]$, we first rewrite Eq.~(\ref{243}) as
\begin{equation} \label{stiffnew2}
[X]^T[\delta K][X_a]=[U]-[X]^T[K_o][X_a]
\end{equation}
As before, Eq.~(\ref{stiffnew2}) can be manipulated
into the following form
\begin{equation} \label{stiffnew4}
[\,B\,] \, \delta\bf{k}=\bf{h}
\end{equation}
where 
\begin{equation} \label{dk}
\delta{\bf k}=\left[
\delta k_{11} \cdots \delta k_{1N} \;|\;
\delta k_{21} \cdots \delta k_{2N} \;|\; \cdots \;|\;
\delta k_{N1} \cdots
\delta k_{NN} \right]^T
\end{equation}
and 
\begin{equation} 
{\bf h}=\left[
h_{11} \cdots h_{1N_e} \;|\;
h_{21} \cdots h_{2N_e} \;|\; \cdots \;|\;
h_{N_e1} \cdots
h_{N_eN_e} \right]^T
\end{equation}
In Eq.~(\ref{dk}), $\delta k_{ij}$ represents
the $(i,j)$th element of $[\delta K]$.
Matrix $[B]$ is of size $N_e^2\times N^2$
and vector $\bf h$ is of length
$N_e^2$.  The components of 
$[B]$ and $\bf h$ can be obtained by expanding the
left and right hand sides of Eq.~(\ref{stiffnew2}), respectively.



\subsection{Solution Technique} \label{solutiontechnique}

%LdeP NOTE
%NOTE -- [[THEOREM A]]: HERE WE ADD A SIMPLE THEOREM PROVING THAT $[G]$ in fact
%has full rank (whether or not the system is overdetermined or underdetermined
%or square).  Then we can safely claim that the overdetermined system has
%one solution, whereas the underdetermined system as infinite solutions.

Equations~(\ref{massnew}) and (\ref{stiffnew4}) are of the general form
\begin{equation} \label{both}
[\,G\,]\,{\bf y}={\bf z}
\end{equation}
where matrix $[\,G\,]$ and vector ${\bf z}$ are both
known and of size $N_e^2\times N^2$ and length $N_e^2$, respectively.
Because the modal matrices $[X]$ and $[X_a]$ % and $[X_b]$
are always of full rank, matrix $[G]$ will also 
be of full rank (see Appendix A for 
a detailed proof).  
When $N_e=N$, Eq.~(\ref{both}) can be solved exactly by 
using simple Gaussian elimination.
When $N_e<N$, Eq.~(\ref{both}) yields an underdetermined
problem (that is, the number of
equations is less than the number of unknowns), 
which, because the system has full rank, will have
an infinite number of solutions. % [reference].  
To render the solution unique,
we may choose a solution vector ${\bf y}$
such that the Euclidean norm of vector ${\bf y}$
is minimized.  The resulting solution is referred to as the
minimum norm least squares solution to Eq.~(\ref{both}).  
%If
%the least squares problem has more than one solution, as in the case
%of an underdetermined system, the
%solution vector having the minimum Euclidean norm is called the minimum norm
%solution.  
Because the
analytical and the actual system matrices are presumed to be close,
%then if there is an infinite set of solutions $\delta{\bf m}$ and $\delta{\bf k}$
%to equations (\ref{massnew}) and (\ref{stiffnew4}),
it is reasonable to use the unique minimum norm solutions 
%$\delta{\bf m}$ and $\delta{\bf k}$
to update the analytical mass and stiffness matrices, respectively.

Initially it may appear that one must solve two
underdetermined, least squares problems of size $N_e^2\times N^2$
(assuming $N_e<N$)
in order to update the system matrices [see Eqs.~(\ref{massnew}) and
(\ref{stiffnew4})].  However,
the optimal matrix storage scheme
commonly used in finite elements \cite{Bathe} can be applied
to pass along the available sparsity information of the
analytical system, and to impose
the condition that all zero elements in the analytical 
system matrices remain zeros in the adjusted system 
matrices.  
%That is, $[\delta M]$ and
%$[\delta K]$ maintain the same sparsity pattern as
%$[M_o]$ and $[K_o]$, respectively.  
Mathematically, this can be
achieved by eliminating all of the known zero elements from ${\bf y}$
and by deleting all the corresponding columns in $[\,G\,]$.
This drastically reduces the size
of the problem to be solved.  

To illustrate how the connectivity information can
be used to reduce the size of the least squares problem, 
consider a system whose analytical mass matrix is 
diagonal.  Then $\delta m_{ij}=0$ for $i\ne j$, and  
Eq.~(\ref{massnew}) reduces to
\begin{equation} \label{mp}
[\,A'\,] \, \delta{\bf m'}={\bf r}
\end{equation}
where $[\,A'\,]$ is obtained from $[\,A\,]$ by deleting all the
columns that multiply by $\delta m_{ij}$ for $i\ne j$, and
\begin{equation} 
\delta{\bf m'}=\left[
\delta m_{11} \quad \delta m_{22} \;\; \cdots \;\; \delta m_{NN} \right]^T
\end{equation}
Thus, the initial problem of size $N_e^2\times N^2$ 
is reduced to one of size $N_e^2\times N$. 
The resulting least squares problem will be either overdetermined 
(that is, the number of equations is greater than or equal to the 
number of unknowns) or 
underdetermined, depending on whether $N_e^2\ge N$ or
$N_e^2<N$, respectively.  


Similarly, the connectivity information of the 
analytical stiffness matrix can also be enforced to 
reduce the size of the least squares
problem to be solved.  
For instance, if the analytical stiffness matrix is 
tri-diagonal, then $\delta k_{ij}=0$
for $|i-j|>1$, and  
Eq.~(\ref{stiffnew4}) reduces to
\begin{equation} \label{sp}
[\,B'\,] \, \delta{\bf k'}={\bf h}
\end{equation}
where $[\,B'\,]$ is obtained from $[\,B\,]$ by deleting all the
columns that multiply by $\delta k_{ij}$
for $|i-j|>1$, and 
\begin{equation} 
\delta{\bf k'}=\left[
\delta k_{11} \quad \delta k_{12} \;|\; 
\delta k_{21} \quad \delta k_{22} \quad \delta k_{23} \;|\;
\cdots \;|\; \delta k_{NN-1} \quad \delta k_{NN} \right]^T
\end{equation}
Thus, the initial problem of size
$N_e^2\times N^2$ is reduced to one of size $N_e^2\times (3N-2)$.


%We re-emphasize that imposing the connectivity information on the mass
%and stiffness matrices of the system not only preserves
%the physical configuration of the structure, it also 
%substantially reduces the size of the least squares problems
%to be solved. 



The proposed updating algorithms rely on the correctness of the
connectivity information. Because the basis of model updating 
is the analytical model,
the analytical model must capture certain physical attributes
of the actual system.  
%Here, we assume that the analytical model
%and the actual system share the same sparsity pattern.  
%The structural mass and stiffness parameters of the
%actual system, however, may deviate substantially
%from those of the analytical model.
%
%
%
%One advantage of the proposed updating algorithm is that it
%is very forgiving when it comes to deviations of the analytical
%system from the physical structure.  That is, the particular mass
%and stiffness parameters may vary substantially between the two
%systems. 
%
%
Here, we assume that the connectivity information
in the analytical model is correct, and the proposed 
updating algorithms
maintain that physical sparsity pattern in the mass
and stiffness matrices.
It will be shown that the proposed updating schemes
are very forgiving when the mass and stiffness
parameters vary substantially between the analytical and the
actual systems. 
%it comes to deviations of the analytical
%system from the physical structure.  That is, the particular mass
%and stiffness parameters may vary substantially between the two
%systems.  





%&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
%Clearly the analytical model for a given
%structure can vary from one analysis to
%another and will vary with the number and placement of the degrees-of-freedom
%used.  Thus, starting with different models for the same
%structure will result in different final corrected models.  Nevertheless,
%for a given analytical model, the proposed updating algorithms will
%return a corrected model that is more accurate than it was initially.
%
%&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&









 

\section{Numerical Issues} %LdeP Added 4/5/99

%{\bf Lisette:}  Consider the following
%\begin{itemize}
%\item  compare underdetermined and overdetermined problems; numerical
%issues and geometric interpretations
%\item  uniqueness and least squares solutions
%\item  least squares, minimum norm solutions for underdetermined problems
%\item  use the inequalities that we had
%obtained earlier for when the system becomes underdetermined
%and when it becomes overdetermined, give
%a criterion on the minimum number of measured modes
%needed in order to make the problem overdetermined, 
%and hence the solution
%unique.  For instance, in updating the mass matrix,
%we may need $N_e>\sqrt{N}$ measured modes;
%in updating the stiffness matrix, we may need
%$N_e>\sqrt{3 N - 2}$ measured modes (assuming that the
%stiffness matrix is tri-diagonal)
%\item  sensitivity of the least squares problems to perturbations
%\item  accuracy versus the number of measured modes
%\end{itemize}

The numerical issues encountered when solving a 
least squares problem differ in
nature depending on whether the problem is 
overdetermined or underdetermined.
Because matrix $[G]$ of Eq.~(\ref{both}) is of full rank,
a least squares solution to Eq.~(\ref{both}) always exists. % [reference]. 
However, since
the reduced systems of Eqs.~(\ref{mp}) and (\ref{sp}) 
were obtained by deleting certain columns of $[G]$ to enforce the
sparsity information, matrix $[G]$ may suffer rank-deficiency.   
While there exist
techniques for determining the numerical rank of such
systems and subsequently finding the unique minimum norm solution,
numerical perturbations in the data, either in the form of measurement
or round-off errors, can 
still give rise to computationally
induced inaccuracies. 
%
%Let $\hat{\bf x}$ denote the computed least squares solution.  The
%quality of $\hat{\bf x}$ depends on how close $\hat{\bf x}$ is 
%to ${\bf x}$ in the sense of the Euclidean norm measure.
%  
It is therefore necessary to understand how ${\bf y}$ is affected
by perturbations in $[G]$ and ${\bf z}$.

To this end, let us examine the least squares solution of a general
$m\times n$ linear system
%Let ${\bf x}$ denote the
%least squares solution of the $m\times n$ system
\begin{equation} \label{o1}
[A]{\bf x}={\bf b}
\end{equation}
Let $\hat{\bf x}$ represent the
least squares solution of the $m\times n$ perturbed system
\begin{equation} 
[\tilde A] \hat{\bf x} = \tilde{\bf b}
\end{equation}
where 
\begin{equation} 
[\tilde A]=[A]+[\delta A]
\end{equation}
and 
\begin{equation} 
\tilde{\bf b}={\bf b}+\delta{\bf b}
\end{equation}
The perturbations $[\delta A]$ and $\delta{\bf b}$ can reflect
either numerical inaccuracies in the floating point representations 
of $[A]$ and ${\bf b}$, or measurement errors during testing.



A problem is said to be numerically 
{\it ill-conditioned} when small errors
in data lead to a relatively large error in the solution,
regardless of how stable the algorithm is for solving the problem.  
The {\it condition number} of a linear system
reveals the accuracy or inaccuracy of the computed result 
due to small perturbations.  For the linear system of Eq.~(\ref{o1}), 
%
%(where $[A]$ is a rectangular matrix of full rank),
% 
its condition number is given by
\begin{equation} \label{cond}
\cond\left( [A] \right) = ||[A]|| \cdot ||[A]^\dagger||
\end{equation}
where $||[A]||$ denotes the Euclidean norm
of $[A]$ and $||[A]^\dagger||$ is the pseudoinverse
of $[A]$, defined as
%
%\b
%[A]^\dagger = \left([A]^T [A]\right)^{-1} [A]^T \ \ \ 
%\mbox{for} \ \ \ m \geq n
%\end{equation}
%and
%\b
%[A]^\dagger = [A]^T \left([A] [A]^T\right)^{-1} \ \ 
%\mbox{for} \ \ \ m < n
%\end{equation}
%
\begin{equation} 
[A]^\dagger = \left\{ \begin{array}{ll}
\left([A]^T [A]\right)^{-1} [A]^T \;, &  \mbox{for} \; m \geq n
\\ \nonumber
[A]^T \left([A] [A]^T\right)^{-1} \;, &  \mbox{for} \; m < n
\end{array}
\right.
\end{equation}













\subsection{Overdetermined Systems}\label{Overdetermined}

If a general $m \times n$ linear system of Eq.~(\ref{o1}) 
is overdetermined ($m \ge n$) and full rank, 
the least squares solution to the system is {\it unique}.  
Premultiplying Eq.~(\ref{o1}) by $[A]^T$ yields
\begin{equation} 
[A]^T [A] {\bf x} = [A]^T {\bf b}
\end{equation}
Because $[A]$ is full rank and $m \ge n$, $([A]^T[A])^{-1}$ exists. % [reference].
Thus, the unique least squares solution of Eq.~(\ref{o1}) is given by
\begin{equation}  
{\bf x} = \left([A]^T [A]\right)^{-1} [A]^T {\bf b} \label{OLSsol}
\end{equation}
where ${\bf x}$ is the best solution 
in the sense of minimizing $||[A]{\bf x}-{\bf b}||$.


%
%\subsubsection{Perturbing ${\bf b}$}
%
%
%Assuming $[A]$ is exact and perturbing only the right hand side of the
%system from ${\bf b}$ to $\tilde{\bf b}$, where
%$\tilde{\bf b} = {\bf b} + \delta{\bf b}$, then Eq.~(\ref{o1})
%becomes
%\b
%[A] \hat{\bf x} = \tilde{\bf b}
%\e
%where $\hat{\bf x}$ denote the computed least squares solution
%to $\tilde{\bf b}$.  The
%quality of $\hat{\bf x}$ depends on how close $\hat{\bf x}$ is 
%to ${\bf x}$ in the sense of the Euclidean norm measure.
%Let ${\bf b}_R$ and $\delta {\bf b}_R$ be the respective projections
%of ${\bf b}$ and $\delta{\bf b}$ onto $R\,([A])$.  The 
%Right Perturbation Theorem (REFERENCE: Datta p. 324) states that
%as long as $||{\bf b}_R|| \ne 0$, then
%\b
%\frac{||\hat{\bf x} - {\bf x}||}{||{\bf x}||} 
%\leq \cond\left( [A] \right) \frac{||\delta {\bf b}_R ||}{|| {\bf b}_R ||}
%\e
%where ${\bf x}$ and $\hat{\bf x}$ are the unique solutions of 
%$[A]{\bf x}={\bf b}$ and $[A]\hat{\bf x} = \tilde{\bf b}$. 
%The implication here is that if $\cond\left( [A] \right)$ is
%large, then even small perturbations in $\delta{\bf b}$ may lead
%to significant inaccuracies in the computed least squares solution.
%
%\subsubsection{Perturbing $[A]$}
%
%
%Suppose now that ${\bf b}$ is exact, but that matrix $[A]$
%has been perturbed by $[\delta A]$.  Then Eq.~(\ref{o1}) becomes
%\b
%[\tilde A] \hat{\bf x} = {\bf b}
%\e
%where $[\tilde A] = [A] + [\delta A]$. 
%Let $[\delta A_R]$ and $[\delta A_o]$ represent the projection of
%$[\delta A]$ onto $R(A)$ and 
%the projection of $[\delta A]$ 
%onto the orthogonal complement of $R\,([A])$, respectively.
%Then Stewart (REFERENCE: see Datta p.326) has shown that if 
%Let ${\bf x}$ and $\hat{\bf x}$ be the unique solutions of 
%$[A]{\bf x}={\bf b}$ and $[\tilde A]\hat{\bf x} = {\bf b}$, 
%and assume that $[\delta A]$ is small enough so that
%\b
%\mbox{rank}\,([A]) = \mbox{rank}\,([A] + [\delta A])
%\e
%then if $|| {\bf b}_R ||\ne 0$ [reference]
%\b\label{46}
%\frac{||\hat{\bf x} - {\bf x}||}{||{\bf x}||} \leq 2\ \cond\left( [A] \right) 
%\frac{|| [\delta A_R] ||}{|| [A] ||}
%+ 4\ \left( \cond [A] \right)^2 
%\frac{|| [\delta A_o]] ||}{||[A] ||} \frac{|| {\bf b}_o ||}{|| {\bf b}_R ||}
%    + {\cal O}\left(\frac{|| [\delta A_o] ||}{||[A] ||}\right)^2
%\e
%where ${\bf b}_o$ is the projection of ${\bf b}$ onto
%the orthogonal complement of $R\,([A])$.
%The implication of Eq.~(\ref{46}) is that as long as $|| [\delta A_o] ||$
%and $|| {\bf b} ||$ are nonzero, then the accuracy
%of the least squares solution depends in general on the {\it square} of 
%$\cond\left( [A] \right)$.  If either $|| [\delta A_o] ||$
%or $|| {\bf b} ||$ is zero, the numerical sensitivity of the solution
%depends only on $\cond\left( [A] \right)$.
%
%
%\subsubsection{Perturbing $[A]$ and ${\bf b}$}
%

The effects of perturbing $[A]$ and ${\bf b}$ separately
on the solution of Eq.~(\ref{o1}) have been well studied and
are detailed in \cite{Datta}.
%on the solution of Eq.~(\ref{o1}) are detailed, 
%for example, Datta \cite{datta}.
The results become more complicated when 
$[A]$ and ${\bf b}$ are perturbed simultaneously.  
Golub and van Loan in \cite{GvL} %\cite[p.228]{GvL} %(REFERENCE: GVL p 228) 
provide the following theorem that can be used to analyze
the condition of the system when both
$[A]$ and ${\bf b}$ undergo perturbations.
%
%Suppose ${\bf x}$ and $\hat{\bf x}$ are the respective 
%least squares solutions of the
%$m\times n$ ($m \geq n$) systems $[A]{\bf x}={\bf b}$ 
%and $[\tilde A] \hat{\bf x} = \tilde{\bf b}$,
%and $|| {\bf b} || \neq 0$.  
%
Defining
\begin{equation} \label{ep}
\epsilon = \max\left\{\frac{|| [\delta A] ||}{|| [A] ||},
\frac{|| \delta{\bf b} ||}{|| {\bf b} ||}   \right\}
\end{equation}
and letting $\theta$ be such that (for $|| {\bf b} || \neq 0$)
\begin{equation}  
\sin{(\theta)} = \frac{|| [A]{\bf x} - {\bf b} ||}{|| {\bf b} ||} \neq 1
\end{equation}
then the upper bound on the relative solution error is given by
\begin{equation} 
\frac{||\hat{\bf x} - {\bf x} ||}{|| {\bf x} ||} 
\leq \epsilon \left( \frac{2\ \cond([A])}{\cos{(\theta)}} 
+ \tan{(\theta)} \cdot \left( \cond([A]) \right)^2 \right) 
+ {\cal O}(\epsilon^2) \label{OLSsolerror}
\end{equation}
where
\begin{equation}  
\tan{(\theta)} = \frac{|| [A]{\bf x}-{\bf b} ||}
{\sqrt{|| {\bf b} ||^2-|| [A]{\bf x}-{\bf b} ||^2}}
\end{equation}
Thus, if the residual norm $|| [A]{\bf x}-{\bf b}||$ is nonzero, 
it is the {\it square} of the $\cond([A])$ which measures 
the sensitivity of the least squares solution
to numerical perturbations.

%
%\subsubsection{Numerical Solution Approaches for Overdetermined Systems}
%

The Normal Equations Method %[cf, REFERENCE] 
and the QR factorization algorithms %[cf, REFERENCE] 
are commonly used to solve an overdetermined least 
squares problem (for example, see \cite{Datta}).  
Of the QR methods, Householder transformations, Givens rotations,
or modified Gram-Schmidt orthogonalizations %[references] 
can be implemented.
QR factorizations tend to be more numerically reliable, but are more expensive,
in terms of total FLOP counts, by at least a factor of two. 
According to Golub and van Loan \cite{GvL}, %\cite[p.230]{GvL}, %(REFERENCE: GvL p 230), 
the Normal Equations solution
method produces an $\hat{\bf x}$ whose relative error always depends on the {\it square}
of $\cond([A])$, whereas the QR approaches produce solutions which depend on both
$\cond([A])$ and the product $|| [A]{\bf x}-{\bf b}|| \cdot ( \cond[A])^2$.  
Thus, if the norm of 
the residual $|| [A]{\bf x}-{\bf b}||$ is small, 
QR has a distinct advantage over Normal Equations. 
However, if $|| [A]{\bf x}-{\bf b}||$ is large 
and the problem is ill-conditioned, both methods
will produce inaccurate results.
Finally, because the extra computational expense necessary for implementing
the QR solution method is a relatively small 
price to pay for increased numerical stability, 
a QR approach is used in this paper to solve overdetermined systems.



\subsection{Underdetermined Systems}\label{Underdetermined}


The fundamental difference between the solution 
of an overdetermined and an underdetermined
least squares problem of the form of Eq.~(\ref{o1})
of full rank is that while there is a {\it unique} least
squares solution to the overdetermined system ($m \ge n$), 
there are an {\it infinite} number
of solutions to the underdetermined system ($m < n$).% [reference].
  

This can be understood when we consider
that if $[A]{\bf x}={\bf b}$ is an underdetermined system of equations,
%(that is, of dimension
%$m\times n$ with $m < n$), 
then every least squares solution $x$ can be represented
as
\begin{equation} 
{\bf x} = {\bf x}_p + {\bf x}_N
\end{equation}
where ${\bf x}_p$ is any solution satisfying $[A] {\bf x}_p = {\bf b}$, 
and ${\bf x}_N$ is in the 
null space of $[A]$.  That is,
\begin{equation} 
[A] {\bf x}_N = {\bf 0}
\end{equation}
The complete solution set is given by %\cite[p.349]{Datta}%(REFERENCE: cf, Datta p 349)
\begin{equation} 
{\bf x} = [A]^\dagger {\bf b} + \left( [I] - [A]^\dagger [A]\right) {\bf y} \label{ULSsol}
\end{equation}
where ${\bf y}$ is an arbitrary $n$-vector.  Recall that since $m < n$, 
$[A]^\dagger$ is defined by $[A]^\dagger = [A]^T ([A] [A]^T)^{-1}$.
Comparing Eqs.~(\ref{OLSsol}) and
(\ref{ULSsol}), we note the absence of the
arbitrary vector, ${\bf y}$, in Eq.~(\ref{OLSsol}).

%
%
%\subsubsection{Perturbing $[A]$ and ${\bf b}$}
%

When there are an infinite number of solutions, it is commonly the
minimum norm solution that is sought and then used to correct
the system matrices.  Additionally,
perturbation analysis has been done mainly for the minimum norm case.
%
%The numerical sensitivity results for an underdetermined
%problem are found to be similar to those for an overdetermined problem.  
%Suppose ${\bf x}$ and $\hat{\bf x}$ 
%are the respective {\it minimum norm} least squares solutions 
%of the
%$m\times n$ ($m \leq n$) systems $[A]{\bf x}={\bf b}$ 
%and $[\tilde A] \hat{\bf x} = \tilde{\bf b}$,
%and ${\bf b} \neq {\bf 0}$.  Let
%\b
%\epsilon = \max\left\{\frac{|| [\delta A] ||}{|| [A] ||},
%\frac{|| \delta{\bf b} ||}{|| {\bf b} ||}   \right\}
%\end{equation}
%
Golub and van Loan \cite{GvL}
%(REFERENCE: GvL p.258) 
show that for an underdetermined
system of full rank, an upper
bound on the relative error in the minimum norm solution is given by
\begin{equation} 
\frac{|| \hat{\bf x} - {\bf x} ||}{|| {\bf x} ||} \leq \cond( [A] ) \ 
\left(\frac{|| [\delta A] ||}{|| [A] ||}\min{(2,n-m+1)} + 
\frac{|| \delta{\bf b} ||}{|| {\bf b} ||} \right)
  + {\cal O}(\epsilon^2) \label{ULSsolerror} 
\end{equation} 
%
%where $\cond([A])$ is still given by Eq.~(\ref{cond}), but
%now the pseudoinverse of $[A]$ is defined as
%\b
%[A]^\dagger = [A]^T \left([A] [A]^T\right)^{-1} \ \ \mbox{for} \ \ \ m < n
%\end{equation}
%
Equation (\ref{ULSsolerror}) shows that
the sensitivity of an underdetermined system now depends
on $\cond( [A] )$, as opposed to the overdetermined case in which
the sensitivity depends on $(\cond [A] )^2$ (see Eq.~(\ref{OLSsolerror})).
Because $\cond([A])\geq 1$ for any matrix $[A]$, the
error bound for the overdetermined case can become fairly large unless
the norm of the residual is sufficiently small.  However, in the underdetermined
case, the error bound does not depend explicitly on the norm of the residual,
but only on the perturbations in $[A]$ and $\bf b$.

%
%\subsubsection{Numerical Solution Approaches for Underdetermined Systems}
%

Algorithms for finding the least squares solution of an
underdetermined system include the
Normal Equations method and the QR factorization algorithms.  
A solution approach
may either find some non-unique least squares solution to $[A]{\bf x}={\bf b}$, or
find the unique solution ${\bf x}$ such 
that the Euclidean norm of ${\bf x}$, $|| {\bf x} ||$, is minimized.
The latter, which is the minimum norm solution, is given explicitly
by 
\begin{equation} 
{\bf x} = [A]^\dagger {\bf b}
\end{equation}
Of course, it is entirely possible that 
another solution besides the one with minimum norm 
may yield a more accurate updated model.   
However, as discussed previously,  because $[\delta M]$ and $[\delta K]$ 
are assumed to be relatively small, it
is reasonable to use the unique minimum norm solution in this case
as corrections to the initial analytical matrices.
%In this case, however, it is not
%practical to search the entire solution set for the best fit, so we solve
%for the unique solution with minimum norm.

%that the updated model obtained by solving
%an underdetermined system may not be as accurate as that
%obtained by solving an overdetermined system.  



\subsection{Relating Accuracy to $N_e$ and $N$}

For systems in which there is no special sparsity pattern to be maintained,
the least squares system will always have dimension $N_e^2 \times N^2$
[see Eqs.~(\ref{massnew}) and (\ref{stiffnew4})].
If the measured modes are incomplete (for $N_e<N$), the problem
will always be underdetermined.  
%In this case we note that the closer
%$N_e$ is to $N$, the more accurate our solutions are likely to be.  
%Replacing $m$ by $N_e^2$ and $n$ by $N^2$ in the perturbation
%bound given in Eq.~(\ref{ULSsolerror}), we observe that the closer $N_e$ is to
%$N$, the smaller the factor on $\cond( [A] )$ is ({\bf not
%quite true: since we're taking $\min{(2,n-m+1)}$, for
%$n-m+1>2$, the factor will always be 2!}), and thus the smaller the bound
%on our minimum norm least squares solution error.  


As stated previously, an overdetermined system of full rank
has a {\it unique} least squares solution, whereas an 
underdetermined system of full rank admits a unique solution
only when the additional minimum norm constraint is imposed. 
Often the physical structure will possess a special
sparsity pattern in the analytical system that we will have to maintain.
%When the sparsity information is imposed, the system may suffer rank deficiency.  
Our numerical experiments indicate 
that in certain cases, when the sparsity pattern
of a physical structure is imposed, 
more accurate solutions are achieved when the least squares system
is overdetermined as opposed to underdetermined, even
if the systems are rank-deficient.  
This may be surprising in light of the fact that our error bounds
(\ref{OLSsolerror}) and (\ref{ULSsolerror}) seem to indicate that there is more room for
error in the overdetermined case.  However, these are upper bounds and
not exact error measurements.  Additionally, if the norm of the residual vector
is small compared to the perturbations in $[A]$ and $\bf b$,
the error bound for the overdetermined system could in fact be smaller than
that of the underdetermined system.
Also, these experimental results make
some intuitive sense: the more information we can gather about the physical system,
the better our updated model becomes. 
Thus, to improve the accuracy of the updated model, 
solving an overdetermined system for this particular application is preferred.

The above discussion has direct implications about
the minimum number of experimentally determined modes, $N_e$, 
which should be measured before implementing the proposed model
updating algorithm.
%
%Suppose the sparsity pattern of the structure is such that the dimension of our
%least squares system becomes $m\times n$.  To induce 
%an overdetermined system requires
%\b
%m \ge n
%\end{equation}
%
For example, if $\delta M_{ij}=0$ for $i\neq j$ 
(that is, $[M]$ is a diagonal matrix),
then the least squares system has dimension $N_e^2 \times N$.  In this
case, to induce an overdetermined system requires a minimum of
\begin{equation} 
N_e \ge \sqrt{N}
\end{equation}
measured modes.  Similarly, in the case of the 
stiffness matrix, if the condition
that $[K]$ is a tridiagonal matrix is imposed, at least
\begin{equation} 
N_e \ge \sqrt{3N-2}
\end{equation}
experimental modes are needed
to achieve an overdetermined least squares system.
The above results will be validated 
the following sections.


\section{Numerical Experiments} \label{numexps}

To update the analytical mass and stiffness matrices using the
proposed routines requires the solution of least squares
problems, which can be either overdetermined or underdetermined,
depending on the number of measured modes, $N_e$, used to
perform the update.  When solving a least squares problem, the CMLIB
\cite{gams} routine
{\it sglss} was implemented, which is specialized to handle 
both underdetermined and overdetermined systems of the form
$[A]{\bf x} = {\bf b}$, where $[A]$ is an $m\times n$ matrix
and ${\bf b}$ is a vector of length $m$.
When the system is overdetermined ($m\ge n$), 
the least squares solution is computed by
decomposing the matrix $[A]$ into the product of an
orthogonal matrix $[Q]$ and an upper triangular matrix
$[R]$ (QR factorization).
When the system is underdetermined $(m<n)$, the minimum
norm solution is computed by factoring the
matrix $[A]$ into the product of a lower triangular
matrix $[L]$ and an orthogonal matrix $[Q]$ (LQ factorization). 
If matrix $[A]$ is determined to be rank
deficient, that is, if the rank of $[A]$ is less than
min($m$,$n$), then again the minimum norm least squares
solution is computed. 




%\subsection{The Effects of Perturbation 
%on the Stability of the Updating Algorithms} \label{stability}
%\begin{itemize}
%\item  In general, the perturbations can arise due to
%either numerical inaccuracies in the floating point
%representations or measurement errors during testing.
%Since the updated model is only as accurate as the measured
%data, it will be assumed that the measured modes
%of vibration are exact.  Thus in this paper,
%the perturbations will be the results of round-off errors only.
%
%
%
%\item  Here, the effects of round-off errors on the robustness
%of the updating algorithms are investigated.  While the
%measured data are assumed to be correct, due to the physical
%limitations of the instrumentations, only finite precision
%can be obtained, {\it i.e.}, experimental measurements
%can only be made to a given number of significant digits, denoted
%by $s$.  Incidentally, all the previous results
%were obtained by using $s=6$, {\it i.e.}, assuming
%all the measurements can be obtained to the 6th significant digit.
%
%
%
%\item  Figure 3 shows the variations of
%$\ep_m$, $\ep_k$ and $\ep_{\omega}$
%as a function of $s$, for $N_e=N$.
%The system parameters are identical to those of Tables 1 and 2.
%Note that $\ep_m$, $\ep_k$ and $\ep_{\omega}$ decrease as 
%$s$ increases.  Nevertheless, even for $s=3$, the updated
%model is considerably more accurate than the initial
%analytical model.
%
%
%
%\end{itemize}

\subsection{Model Updating Based on a 
Complete Set of Measured Modes} \label{complete}


Consider the system of Fig.~\ref{Figure1} with $N=25$.  
We first examine the effectiveness of
our model updating techniques for the case where the number of measured modes,
$N_e$, is equal to the total degrees-of-freedom of the system, $N$ (that is,
the set of measured modes is complete).  The analytical 
masses and stiffnesses are given by $m_o = 2.00$ kg and
$k_o = 5.00$ N/m.  The actual masses and stiffnesses
are listed in Tables 1 and 2, respectively.

%\begin{figure} 
%{\epsfig{figure=oscillators.eps}}
%\caption{Simple chain of coupled oscillators.}
%\label{Figure1}
%\end{figure}
\begin{figure} 
\centering\scalebox{0.6}[0.9]{\includegraphics{oscillators.eps}}
\caption{Simple chain of coupled oscillators.}\label{Figure1}
\end{figure}
%%\centerline{Figure 1:  Simple chain of coupled oscillators.}


\tabcap{2in}
\begin{table} \label{Table1}
\caption{The actual and the
updated masses (in kg), for $N_e=N=25$.
The analytical masses are $m_o=2.000$ kg.}
\begin{center}
\begin{tabular}{|c|c||c|c|} \hline
$m_{\rm actual}$ &  $m_{\rm update}$ & $m_{\rm actual}$ & $m_{\rm update}$  \\ \hline 
$m_1=$1.2942 &   1.2942 &  $m_{14}=$2.6722 &   2.6722 \\ \hline
$m_2=$2.1831 &   2.1831 &  $m_{15}=$1.3355 &   1.3355 \\ \hline
$m_3=$1.3117 &   1.3117 &  $m_{16}=$2.6680 &   2.6680 \\ \hline
$m_4=$2.6581 &   2.6581 &  $m_{17}=$1.3899 &   1.3899 \\ \hline
$m_5=$2.4371 &   2.4371 &  $m_{18}=$1.8632 &   1.8632 \\ \hline
$m_6=$1.7651 &   1.7651 &  $m_{19}=$1.6231 &   1.6231 \\ \hline
$m_7=$2.8502 &   2.8502 &  $m_{20}=$1.1578 &   1.1578 \\ \hline
$m_8=$1.7984 &   1.7984 &  $m_{21}=$1.1439 &   1.1439 \\ \hline
$m_9=$1.7793 &   1.7792 &  $m_{22}=$1.8189 &   1.8189 \\ \hline
$m_{10}=$2.7588 &   2.7588 &  $m_{23}=$1.2320 &   1.2320 \\ \hline
$m_{11}=$2.1221 &   2.1221 &  $m_{24}=$1.6389 &   1.6389 \\ \hline
$m_{12}=$1.1613 &   1.1613 &  $m_{25}=$2.2254 &   2.2254 \\ \hline
$m_{13}=$2.0234 &   2.0234 &   -----          &   -----  \\ \hline
\end{tabular}
\end{center} 
\end{table}
%\item{Table 1.  The actual and the
%updated masses (in kg), for $N_e=N=25$.
%The analytical masses are $m_o=2.000$ kg.} 




%\newpage



\tabcap{2in}
\begin{table} \label{Table2}
\caption{The actual and the
updated stiffnesses (in N/m), for $N_e=N=25$.
The analytical masses are $k_o=5.000$ N/m.}
\begin{center}
\begin{tabular}{|c|c||c|c|} \hline
$k_{\rm actual}$ &  $k_{\rm update}$ & $k_{\rm actual}$ & $k_{\rm update}$  \\ \hline 
$k_1=$4.1400 &   4.1398 &  $k_{14}=$5.9877 &   5.9877\\ \hline
$k_2=$6.8802 &   6.8804 &  $k_{15}=$5.4973 &   5.4973\\ \hline
$k_3=$5.6052 &   5.6056 &  $k_{16}=$5.9483 &   5.9483\\ \hline
$k_4=$6.5108 &   6.5109 &  $k_{17}=$5.0320 &   5.0320\\ \hline
$k_5=$2.9343 &   2.9343 &  $k_{18}=$6.3608 &   6.3608\\ \hline
$k_6=$7.1326 &   7.1326 &  $k_{19}=$6.2726 &   6.2726\\ \hline
$k_7=$3.3072 &   3.3072 &  $k_{20}=$7.0572 &   7.0572\\ \hline
$k_8=$3.2986 &   3.2986 &  $k_{21}=$6.6026 &   6.6026\\ \hline
$k_9=$6.2021 &   6.2020 &  $k_{22}=$4.9326 &   4.9326 \\ \hline
$k_{10}=$6.6399 &   6.6398 & $k_{23}=$6.3932 &   6.3932 \\ \hline
$k_{11}=$5.9489 &   5.9488 & $k_{24}=$5.8004 &   5.8003 \\ \hline
$k_{12}=$6.3203 &   6.3202 & $k_{25}=$6.5935 &   6.5935 \\ \hline
$k_{13}=$3.3357 &   3.3356 & -----           &   -----   \\ \hline
\end{tabular}
\end{center} 
\end{table}



%\newpage

\tabcap{2in}
\begin{table} \label{Table3}
\caption{The analytical, actual and
updated eigenvalues (in 1/sec$^2$),
for the analytical, actual and updated system
parameters of Tables \ref{Table1} and \ref{Table2}}
\begin{center}
\begin{tabular}{|c|c|c||c|c|c|} \hline
$\lambda_{\rm analytical}$ & $\lambda_{\rm actual}$ 
&  $\lambda_{\rm update}$ & $\lambda_{\rm analytical}$ 
& $\lambda_{\rm actual}$  &  $\lambda_{\rm update}$ \\ \hline 
$\l_1=$0.0095 &0.0105 &   0.0105 &  $\l_{14}=$5.4613 &6.6635 &   6.6635 \\ \hline
$\l_2=$0.0851 &0.0954 &   0.0954 &  $\l_{15}=$6.0697 &7.7364 &   7.7365 \\ \hline
$\l_3=$0.2353 &0.2580 &   0.2580 &  $\l_{16}=$6.6618 &8.8005 &   8.8005 \\ \hline
$\l_4=$0.4577 &0.4674 &   0.4674 &  $\l_{17}=$7.2287 &9.8139 &   9.8139 \\ \hline
$\l_5=$0.7489 &0.8441 &   0.8441 &  $\l_{18}=$7.7618 &9.9112 &   9.9112 \\ \hline
$\l_6=$1.1046 &1.1282 &   1.1282 &  $\l_{19}=$8.2531 &10.7379 &  10.7378 \\ \hline
$\l_7=$1.5193 &1.7577 &   1.7577 &  $\l_{20}=$8.6950 &11.3051 &  11.3051 \\ \hline
$\l_8=$1.9868 &2.1220 &   2.1220 &  $\l_{21}=$9.0810 &12.2057 &  12.2056 \\ \hline
$\l_9=$2.5000 &2.7864 &   2.7864 &  $\l_{22}=$9.4051 &12.9773 &  12.9774 \\ \hline
$\l_{10}=$3.0511 &3.2055 &   3.2056 & $\l_{23}=$9.6624 &13.2552 &  13.2555 \\ \hline
$\l_{11}=$3.6317 &3.5768 &   3.5768 & $\l_{24}=$9.8490 &14.8620 &  14.8619 \\ \hline
$\l_{12}=$4.2330 &4.9515 &   4.9515 & $\l_{25}=$9.9621 &18.8955 &  18.8955 \\ \hline
$\l_{13}=$4.8460 &5.4920 &   5.4920 & -----           &   -----  &   -----  \\ \hline
\end{tabular}
\end{center} 
\end{table}
%\item{Table 3.  The analytical, actual and
%updated eigenvalues (in 1/sec$^2$),
%for the analytical, actual and updated system 
%parameters of Tables 1 and 2.}
%The analytical eigenvalues are given by $(\l_i)_{\rm analytical}=
%4{\displaystyle{\frac{k_o}{m_o}}\sin^2\frac{(2i-1)\pi}{2(2N+1)}  }\;$.} 



For the system of Fig.~\ref{Figure1}, because the mass
matrix $[M]$ is diagonal, updating the mass matrix requires 
finding the solution of
an $N_e^2\times N$ least squares problem.
Because the stiffness matrix is tri-diagonal, updating the stiffness matrix
requires finding the solution of an $N_e^2\times (3N-2)$ 
least squares problem. 
In general, these systems may be overdetermined or
underdetermined, depending on whether 
$N_e^2$ is greater than or less than $N$ in the case of
mass updating, and on whether $N_e^2$ is greater than or
less than $(3N-2)$ in the case of stiffness updating.
For this set of experiments, because $N_e=N$, both the mass update
system and the stiffness update system will be overdetermined.

In Table 1 we compare the actual and the updated masses for
the system of Fig.~\ref{Figure1}.  
For the purpose of numerical simulations, the $i$th actual mass,
$m_i$, and the analytical mass, $m_o$, are related as follows:
\begin{equation} 
m_i=m_o(1+\Delta m_i)
\end{equation}
where $\Delta m_i$ represents the percent deviation of the
$i$th actual mass from its nominal analytical value.  The $\Delta m_i$'s
are randomly chosen (with uniform distribution), and
they have a mean and standard deviation of $-0.06$ 
and 0.28, respectively.
To perform the mass updating
algorithm, lumped masses of magnitude $0.2$ kg are added
to masses (or nodes) 5, 10, 15, 20 and 25.  The added 
masses are all an order of magnitude smaller than the
nominal analytical masses.
%
%, and they
%are chosen so that their magnitude and location satisfy
%the criterion of Eq.~(\ref{a5}).
Note how well the updated masses correspond to the actual values,
despite the large deviations of the actual masses
from the analytical values.


%\newline{\bf Use new stiffness updating algorithm}  

Table 2 shows the actual and the updated stiffnesses
for the system of Fig.~\ref{Figure1}. 
As in the case of the masses, the numerical 
choice for the actual stiffnesses, $k_i$,
and the analytic stiffnesses, $k_o$, are chosen to have the following relationship:
\begin{equation} 
k_i = k_o(1+\Delta k_i)
\end{equation}
where $\Delta k_i$ represents the percent deviation of the
$i$th actual stiffness from its nominal analytical
value.  The $\Delta k_i$'s
are randomly chosen (with uniform distribution), and 
they have a mean and standard deviation of 0.13 and 0.25,
respectively.
%To update the stiffness matrix, grounded springs
%of stiffness $0.2$ N/m are attached 
%to masses (or nodes) 5, 10, 15, 20 and 25.  
%Note that the stiffnesses 
%of the added grounded springs are all an order of magnitude smaller than the
%nominal analytical stiffnesses of the coupling springs.
In this case, too, the updated stiffness values track the 
actual stiffnesses very closely, even though the
deviations between the actual and the analytical stiffnesses
are large.  



In Table 3, we compare the analytical eigenvalues for the system
of Fig.~\ref{Figure1} with those of the actual and the updated systems.  
The analytical eigenvalues of the system
are found by solving Eq.~(\ref{eigan}), the actual eigenvalues
are obtained from a modal survey, and the updated
eigenvalues are calculated by solving 
\begin{equation} 
[K_{\rm update}]{\bf x}_{\rm update}=
\lambda_{\rm update}[M_{\rm update}]{\bf x}_{\rm update}
\end{equation}
where $[M_{\rm update}]$ and $[K_{\rm update}]$ are the updated 
mass and stiffness matrices (see Tables 1 and 2
for the updated system parameters), and 
$(\lambda_{\rm update},{\bf x_{\rm update}})$ are the updated
eigenvalues and eigenvectors.
From Table 3, note the excellent agreement between the updated and the
actual (or the measured) eigenvalues of the system, despite
the large differences between the actual and the initial
analytical eigenvalues.  For $i=1,\cdots, N$, 
$(\l_{\rm actual}-\l_{\rm update})_i$ 
%the absolute error between the actual and the updated system 
is consistently two to four orders of magnitude smaller than the corresponding
$(\l_{\rm actual}-\l_{\rm analytical})_i$, which
clearly indicates the drastic improvement in the eigenvalues
of the updated model over
the initial analytical system. 
%the corresponding absolute error between actual and the analytic system.


%
%\item  In order to maintain numerical stability, we 
%add masses or stiffnesses such that the eigenvalues of
%the modified systems are all different from those
%of the actual system.  This can be achieved by adding
%masses or stiffnesses that are an order of magnitude
%less than the nominal values so that the modified eigenvalues
%are perturbations of the eigenvalues of the actual system.
%\item  The above has beneficial implications in that the added
%lumped elements can be small.  Thus, the updating procedure
%can be can be easily implemented experimentally.
%

When the set of measured modes of 
vibration is complete, {\it i.e.}, when $N_e=N$, 
the updating algorithms can be used to accurately 
correct the analytical mass and stiffness matrices.
Because the proposed updating algorithms allow the
well known and readily available
connectivity information to be enforced, our mass and
stiffness updating schemes generate matrices that not only
produce excellent numerical agreement between the physical and the
updated system, but also preserve the physical configuration of the
structure exactly.

Finally, while 5 lumped masses were
used to generate the results of Tables 1 and 2,
fewer lumped masses can also be used to
execute the updating algorithms.  For example, when only 2 lumped 
masses of magnitude 0.2 kg are added at nodes 1 and 25, the 
resulting updated masses are found to be nearly as accurate as those
of Table 1.  However, %When only 2 grounded
%springs of stiffness 0.2 N/m are attached to nodes 1 and 25, 
the proposed stiffness updating algorithm returns
stiffnesses that are not nearly as accurate as those
of Table 2.  Nevertheless, the
resulting updated stiffnesses are still much closer to the actual
values than the initial analytical ones. 
Clearly the number, placement and magnitude of the added
lumped masses affect the quality of the updates.  
Their influences on the proposed updating algorithms
will be detailed in Section \ref{effects}. %presently. %in a later section. %section 5.4.

\subsection{Model Updating Based on an Incomplete Set 
of Measured Modes} \label{incomplete}


The results of section \ref{complete} %the previous section %in section \ref{complete} 
were obtained for the ideal case
when the set of measured modes is complete.  Note, however,
that because of physical
limitations or time and cost constraints, the set
of measured modes is often incomplete.  Thus in practice,
$N_e$ is almost always less than $N$.  
Of considerable interest,
then, is the effect of the number of measured modes, $N_e$,
on the quality of the updated system matrices.


Theoretically, $N_e$ can include any modes
of vibration of the system.  However, because the lower modes
are typically easier to measure experimentally than the
higher modes, the parameter $N_e$ will be used in the
subsequent analysis to
represent the lowest $N_e$ measured modes of
the structure.  Thus, when $N_e=5$, the first 5 measured modes
will be used to perform the update.



To quantify the accuracy of the mass updating algorithm,
we measure the relative error in the updated masses
with the following relation:
%the following error parameter for the updated masses is introduced:
\begin{equation} 
\epsilon_{m}=\frac{ |{\bf m}_{\rm update} - {\bf m}_{\rm actual}| }
{ |{\bf m}_{\rm actual}| }
\end{equation}
where ${\bf m}_{\rm update}$ and ${\bf m}_{\rm actual}$
are vectors of length $N$ whose elements are the updated
and the actual lumped masses, respectively, and $|{\bf a}|$
represents the Euclidean norm of the vector ${\bf a}$.
To illustrate the improvement of the updated masses
over their initial analytical values, the 
following relative error parameter for the 
analytical masses is also introduced:
\begin{equation} 
(\epsilon_{m})_o=\frac{ |{\bf m}_{\rm analytical} - {\bf m}_{\rm actual}| }
{ |{\bf m}_{\rm actual}| }
\end{equation}
Similar expressions can be defined for the stiffness and 
the eigenvalue error parameters, denoted by 
$\epsilon_k$ and $\epsilon_\l$, respectively.



An error parameter can also be defined for the modal
matrix of the system.  It is common practice to
check for the correctness of the modal matrix by resorting to 
the orthogonality characteristics of the normal modes.  If the
actual modal matrix, $[X]$, is properly normalized, 
then it is orthogonal with respect to the
actual mass matrix, $[M]$, such that
\begin{equation} \label{pp}
[X]^T[M][X]=[I]
\end{equation}
where $[I]$ is the identity.  Because the updated
modal matrix is approximate, replacing the actual
modal matrix, $[X]$, by the updated modal matrix, $[X_{\rm update}]$, 
in Eq.~(\ref{pp}) yields
\begin{equation} \label{triple}
[X_{\rm update}]^T[M][X_{\rm update}]=[I_{\rm update}]
\end{equation}
where $[X_{\rm update}]$ is normalized
so that the diagonal elements of $[I_{\rm update}]$
are identically one.  In general, $[I_{\rm update}]$ is 
a full matrix, but if the update is good, all the off
diagonal terms will be small.
%whose diagonal elements consist of unity.  
%Because the updated
%modal matrix is not exactly orthogonal with respect to $[M]$,
The following can be used to describe the accuracy of the modal matrix 
quantitatively:
\begin{equation} 
\epsilon_{X}=||[I]-[I_{\rm update}]|| 
\end{equation} 
As before, to gauge the improvement of the updated modal matrix
over the initial analytical modal matrix, the following 
error parameter for the analytical modal matrix is also
introduced:
\begin{equation}  
(\epsilon_X)_o=||[I]-[X_{\rm analytic}]^T [M] [X_{\rm analytic}]|| 
\end{equation} 
where $[X_{\rm analytic}]$ 
represents the normalized analytical modal matrix of the
system.


For an updated model to be judged better
than the initial analytical model, we must have
$\epsilon_{m}<(\epsilon_{m})_o$, $\epsilon_{k}<(\epsilon_{k})_o$, 
$\epsilon_{\l}<(\epsilon_{\l})_o$ and $\epsilon_X < (\epsilon_X)_o$.
For an updated model to be considered accurate, 
$\epsilon_{m}$, $\epsilon_{k}$, $\epsilon_{\l}$ and $\epsilon_X$ must be sufficiently small.  
Finally, the smaller the error parameters, the better the updated model.


Figure~\ref{Figure2} shows the variations of $\epsilon_{m}$ and
$\epsilon_{k}$ for the system parameters of Tables 1 and 2 as a function of
$N_e$. 
%the first $N_e$ measured modes.
Also shown are the corresponding $(\epsilon_{m})_o$ and $(\epsilon_{k})_o$,
which are independent of $N_e$
and are given by the horizontal lines.
Note the improvement in the updated parameters as $N_e$ increases.
The experimental results are consistent with
physical intuition:  the larger
the knowledge space or the more information
we can gather about the physical system, the better
our updated model becomes.

\begin{figure} 
\centering\scalebox{0.8}[1.0]{\includegraphics{massstiff5NODES.ps}}
\caption{The mass and stiffness
error parameters, $\epsilon_m$ and $\epsilon_k$, as a function
of $N_e$, for the
system parameters of Tables 1 and 2.  Lumped
elements are added to masses 5, 10, 15, 20 and 25.
The horizontal lines represent the mass and stiffness
error parameters of the analytical model, $(\epsilon_m)_o$
and $(\epsilon_k)_o$.}\label{Figure2}
\end{figure} 



The curves of Fig.~\ref{Figure2} reveal the fewest $N_e$
needed in order to achieve a certain level of accuracy.
Therefore, these curves can be used to determine the smallest $N_e$ 
that should be obtained for the purpose of performing
the update.  Interestingly, each curve reaches a saturation point
beyond which additional information does not lead
to significant improvement in the corrected model.  The results clearly
indicate that there is a minimum $N_e$ that
needs to be obtained in order to ensure sufficient accuracy in
the adjusted model.

Figure~\ref{Figure3} shows the resulting error parameters for
the updated eigenvalues and modal matrix,
$\epsilon_{\l}$ and $\epsilon_X$, as a function of $N_e$. 
Note how well the modes of vibration of the updated
system track the actual system, especially as
$N_e$ becomes large.  Thus, the proposed updating
algorithms return system matrices that yield modes
of vibration that are much closer to the measured
data than those of the initial analytical model.

\begin{figure} 
\centering\scalebox{0.8}[1.0]{\includegraphics{EigvalModal5NODES.ps}}
\caption{The eigenvalue and modal matrix
error parameters, $\epsilon_\lambda$ and $\epsilon_{X}$, as a function
of $N_e$, for the
system parameters of Tables 1 and 2.  Lumped
elements are added to masses 5, 10, 15, 20 and 25.
The horizontal lines represent the mass and stiffness
error parameters of the analytical model, $(\epsilon_\lambda)_o$
and $(\epsilon_{X})_o$.}\label{Figure3}
\end{figure}



\subsection{The Minimum $N_e$ Needed to Ensure 
Sufficient Accuracy} \label{minmodes}

An overdetermined system of full rank has a unique least
squares solution.  An underdetermined system of full rank will admit
a unique least squares solution only
when an additional constraint is imposed, such
as the minimum norm constraint.
In our application,
when the sparsity information of the mass or stiffness
matrix is enforced, the system may become rank-deficient.  From
numerical experiments, it was observed that more accurate solutions
are obtained when the least squares system is overdetermined 
as opposed to underdetermined, even if the systems are rank-deficient.
The numerical results from 
%the previous section have direct
section \ref{incomplete} have direct
implications about the minimum number of experimentally determined modes
one should measure in order to ensure sufficient accuracy of the
updated solution.

For instance, to update a diagonal mass matrix so that
the resulting least squares problem becomes overdetermined, we would need
at least $N_e\ge\sqrt{N}$ measured modes.
Similarly, to update a tri-diagonal stiffness matrix so that
the resulting least squares problem is overdetermined, we would need
at least $N_e\ge\sqrt{3N-2}$ measured modes.  


The above criteria regarding the fewest number
of measured modes needed to perform the update are supported by
various numerical experiments.  For $N=25$, our heuristic
indicates that at least 
5 measured modes are required to sufficiently update the mass matrix, and
at least 9 measured modes are needed to correct the stiffness
matrix.  Naturally, because the number of measured modes
used to perform the update directly affects
the errors, $\epsilon_m$ and $\epsilon_k$, the minimum number
of measured modes that we should use to perform the
update also depends on the size of error
that we are willing to tolerate.
In the case of the experiments of Section \ref{incomplete},
%the previous section, %section \ref{incomplete}, 
the minimum number of measured modes suggested
by the heuristics for the mass update ($N_e=5$)
leads to $\epsilon_m\approx 0.11$ (see Fig.~\ref{Figure2}), and
for the stiffness update ($N_e=9$) results
in $\epsilon_k\approx 0.07$ (see Fig.~\ref{Figure2}), both of which are substantially
lower than their corresponding analytical error parameters.


  
It should be emphasized that the criteria regarding
the fewest $N_e$ needed to perform the update
are formulated empirically.  Thus, depending on the system
parameters, and, as mentioned previously, required
error tolerances, sometimes more and other times fewer measured
modes are required to ensure accuracy of the updated model.
For example, for the mass parameters of Fig.~\ref{Figure2}, 
$7$ measured
modes are sufficient to correct the mass matrix
to achieve an error of less than $0.025$ ($\epsilon_m < 0.025$).
For the stiffness values of Fig.~\ref{Figure2}, 13 measured modes
are sufficient to update the stiffness matrix
to achieve an error of less than $0.025$ ($\epsilon_k <0.025$).



\subsection{The Effects of Locations, 
Magnitudes and the Number of Added Elements} \label{effects}

%The proposed model updating algorithms
%require the modes of vibration of the mass-modified and
%the stiffness-modified systems to be measured.
From numerical experiments, it was observed that the location,
magnitude and the number of added lumped masses may affect the
numerical stability of the updating algorithms. 
The objective of the following discussion is to formulate an
expression that can be used as a measure for determining
whether the resulting mass-modified system
allow for a numerically stable updated solution.
Additionally, we will develop an expression that will provide
some guidance in choosing the required magnitude, placement and the 
number of added lumped masses in order to execute
the updating algorithms.

%The objectives of the following discussion
%are two fold.  First, we aim to formulate an
%expression which can be used to guide us in choosing
%the magnitude and placement of the added lumped elements.
%Second, we plan to show how the magnitude and location
%of the added elements may affect the numerical stability
%of the updating algorithm.

%Such effects
%can be understood by using eigenvalue perturbation theory, assuming
%that the magnitudes of the added elements are small compared to the magnitudes
%of the actual elements.


A frequently encountered scenario in structural dynamics 
is determining the changes in the eigensolution of a 
system after certain modifications are introduced.  
If the changes made are small, 
then the initial modal characteristics 
can be used as a basis from which to extract 
the new eigensolution of the modified 
system without performing a new and possibly costly analysis.  
Using perturbation theory for the eigenvalues
and assuming the unperturbed eigenvalues to be distinct, 
it can be shown (see references \cite{courant} 
and \cite{meiro}) that the $j$th first-order perturbed 
eigenvalue of a slightly modified structural system is given by:
\begin{equation} \label{a1}
\rho_j=\rho_{oj}+{\bf y}_{oj}^T \left([\delta{\cal K}] 
- \rho_{oj} [\delta {\cal M}] \right){\bf y}_{oj}
\end{equation}
where $\rho_j$ represents the $j$th eigenvalue of the
perturbed system, $\rho_{oj}$ is the $j$th unperturbed eigenvalue,
${\bf y}_{oj}$ is the $j$th unperturbed eigenvector,
and $[\delta{\cal M}]$ and $[\delta{\cal K}]$ are 
the first-order perturbation
mass and stiffness matrices, respectively.



Consider now the case of correcting or updating
the mass matrix of the system using the 
mass updating %algorithm developed in a previous section to solve Eq.~(\ref{massnew}).
algorithm developed in Section \ref{massupdating} to solve solve Eq.~(\ref{massnew}).
%algorithm developed in section 3.1,
%The objectives of the following discussion
%are two fold.  First, we aim to formulate an
%expression which can be used to guide us in choosing
%the magnitude and placement of the added lumped elements.
%Second, we plan to show how the magnitude and location
%of the added elements may affect the numerical stability
%of the updating algorithm.  Such effects
%can be understood by using eigenvalue perturbation theory.
This method requires a modal survey of the mass-modified system
to be performed.
Assuming that the added masses are small
compared to the analytical masses of the system,
then the unperturbed and the perturbed 
systems correspond to the actual and the
mass-modified systems, respectively.
Using the nomenclature introduced earlier,
we have for this case $[\delta{\cal M}]=[M_a]$ (matrix of added masses), 
$[\delta{\cal K}]=[0]$, $\rho_{oj}=\l_j$ (the $j$th measured
eigenvalue of the actual system) and ${\bf y}_{oj}={\bf x}_j$
(the $j$th measured eigenvector of the actual system).
Then Eq.~(\ref{a1}) can be rewritten as
\begin{equation} 
\l_{aj}=\l_j-\l_j{\bf x}_j^T [M_a]{\bf x}_j \label{eigmass}
\end{equation}
where $\l_{aj}$ denotes the $j$th eigenvalue of the
perturbed or the mass-modified system.
Equation (10) shows that
the mass updating algorithm is stable as long as
\begin{equation} \label{eigstable}
\l_i-\l_{aj}\ne 0, \;\; {\rm for} \;\; i,j=1,\cdots,N_e
\end{equation}
That is, we require all the measured eigenvalues for the mass-modified
system to be distinct from those of the initial
structure.  Using Eq.~(\ref{eigstable}),
each time we add a different set of masses to the
structure, we have to
remeasure eigenvalues of the new system in order
to check our stability criterion.
%
%Thus, each time we add a different set of masses to the actual
%system, we have to recalculate the eigenvalue information for the
%mass-modified system in order to check our stability criterion
%of Eq.~(\ref{eigstable}).
%
We therefore look for an equivalent equation which directly
relates the actual elements
of the added mass system to the stability of the numerical solution.
Substituting Eq.~(\ref{eigmass}) into Eq.~(\ref{eigstable}) we have
\begin{equation} \label{a2}
\l_i-\l_j + \l_j{\bf x}_j^T [M_a]{\bf x}_j \ne 0, \;\; {\rm for} \;\; i,j=1,\cdots,N_e
\end{equation}
Rearranging Eq.~(\ref{a2}), it becomes
\begin{equation} \label{a3}
{\bf x}_j^T [M_a]{\bf x}_j \ne 1 - \frac{\l_i}{\l_j}, \;\; {\rm for} \;\; i,j=1,\cdots,N_e
\end{equation}
Expanding the left hand side of Eq.~(\ref{a3}) yields
\begin{equation} \label{a4}
\sum_{r=1}^{N}\sum_{s=1}^{N} (M_a)_{rs} (x_j)_r (x_j)_s 
\ne 1 - \frac{\l_i}{\l_j}, \;\; {\rm for} \;\; i,j=1,\cdots,N_e
\end{equation}
where $(M_a)_{rs}$ represents the $(r,s)$th element of the
added mass matrix $[M_a]$,
and $(x_j)_r$ denotes the $r$th element of eigenvector ${\bf x}_j$ 
(the $j$th measured eigenvector of the actual system).
When the left and right hand sides of
Eq.~(\ref{a4}) are nearly identical, the mass updating algorithm
becomes unstable.  
To ensure numerical stability, the location, magnitude and the 
number of added masses must be chosen so that the inequality
of Eq.~(\ref{a4}) is maintained.  Thus, Eq.~(\ref{a4}) allows
one to select the size, location and the number of
added masses that should be employed.
%This is accomplished by selecting appropriate 
%values for the added mass matrix elements $(M_a)_{rs}$.
%
%Note that it is not necessary to recalculate eigenvalue and 
%eigenvector information
%for $[M_a]$ each time a modification is made to $[M_a]$ itself.

Assuming $[M_a]$ to be diagonal, Eq.~(\ref{a4}) simplifies to
\begin{equation} \label{a5}
\sum_{r=1}^{N} (M_a)_{rr}(x_j)^2_r \ne 1 - \frac{\l_i}{\l_j}, 
\;\; {\rm for} \;\; i,j=1,\cdots,N_e
\end{equation}
%While the above equation sums over $N$ terms, 
%it should be noted that some of
%the added masses $(M_a)_{ii}$ can be set equal to zero.
%Thus in application, the number of added
%masses may be substantially smaller than $N$, as long as
%$|\sum_{r=1}^{N} (M_a)_{rr}(x_j)^2_r|$ is greater than
%some fraction of $\left|1 - {\l_i}/{\l_j}\right|$,
%the fraction of which can be determined via numerical simulations
%and will be left as a future research topic [{\bf Lisette:  as I recall,
%you had a good discussion here.  Please include.}]
%
%Note that Eq.~(\ref{a5}) is advantageous in the case one wishes
%to make an apriori determination of the number, placement and size of the
%added masses.  However, the computation of the quantities $\l_i/\l_j$,
%($i,j=1,\cdots,N$), may be impractical, especially for large $N$.
%We note that the problem could be recast as a least-squares problem, making
%the computing more reasonable.  
%
%
%
%
Equation (\ref{a5}) can be used 
to make an a priori determination of the number, 
placement and size of the
added masses.  This can be accomplished
by recasting Eq.~(\ref{a5}) as a least-squares problem
of the form $[A]{\bf x}={\bf b}$, where
$[A]$ and ${\bf b}$ are knowns and of dimension
$N_e^2\times N$ and length $N_e^2$, respectively, and vector $\bf x$ is
the unknown, consisting of $[(M_a)_{rr}]_{r=1,\cdots,N}$.  
Computation of the solution  of the resulting
least squares problem may be very expensive, however,
especially for large $N$ and $N_e$.
For our purposes, because
Eq.~(\ref{eigstable}) is much easier to check once all eigenvalue
measurements have been taken, we choose to make use of Eq.~(\ref{eigstable})
as an a posteriori measure of the goodness of the number, placement 
and size of the
added masses.  Our experiments are intended to highlight the fact that
mass placement, number and size do indeed affect the accuracy of our solution.
In all cases considered, we have ensured that Eq.~(\ref{eigstable}) is
satisfied.

In applying the mass updating scheme, the
mechanism that leads to instability will be more complicated
and more difficult to show
analytically because the perturbed eigenvalues of
the mass-modified system include higher-order terms.  Nevertheless,
the above perturbation analysis readily reveals that the magnitude, location
and the number of added masses do in fact affect the stability of the updating
algorithm.

%\newline{\bf Need to generate new results!}

To illustrate how the added lumped masses affect
the quality of the update, consider the actual
system whose masses and stiffnesses are given by
those of Tables 1 and 2.  Instead of adding 5 lumped
masses to execute the updating algorithms as we
did to obtain the results of Fig.~\ref{Figure2}, consider adding only 2
masses of magnitude 0.2 kg each at nodes 3 and 21.  
Figure~\ref{Figure4} shows the variations of the resulting relative
mass error parameter, $\epsilon_m$, and the resulting relative
stiffness error parameter, $\epsilon_k$, as a function of $N_e$.  Note
that $\epsilon_m$ decreases initially with increasing $N_e$ until
$N_e=4$, but increases to above $(\epsilon_m)_o$ suddenly for $N_e=5$ and $N_e=6$.
For $N_e\ge 7$, $\epsilon_m$ becomes nearly zero.
Similarly, note that $\epsilon_k$ decreases with increasing
$N_e$ until $N_e=9$, increases slightly for
$N_e$ equal to $10$ and $11$, and then drops quickly
down to nearly zero for larger values of $N_e$.
Because the relative solution error
for the least squares problem depends on the
condition number of the linear system (recall the discussions
on overdetermined and underdetermined systems in Sections \ref{Overdetermined} 
and \ref{Underdetermined} ) 
%(see sections 4.1 and 4.2), 
we speculate that
the regions in which $\epsilon_m$ and $\epsilon_k$ are above
$(\epsilon_m)_o$ and $(\epsilon_k)_o$ correspond
to those in which the condition number of the least squares
matrix, cond$([A])$, becomes large.  The sensitivity
of an underdetermined system depends on
cond$([A])$, and the sensitivity
of an overdetermined system is generally dominated by
$({\rm cond}([A]))^2$,  so we introduce an instability
measure which equals the logarithm
of ${\rm cond}([A])$ if the system is underdetermined,
and equals the logarithm of ${\rm cond}([A]^2)$ if 
the system is overdetermined.
%that varies as either cond$([A])$ or $({\rm cond}([A]))^2$,
%depending if the problem is underdetermined or overdetermined.
Figure~\ref{Figure5} shows the variation of the instability number as
a function of $N_e$.  When $N=25$, the mass updating algorithm
is underdetermined for $N_e<5$, and overdetermined for $N_e\ge 5$;
similarly, the stiffness updating algorithm
is underdetermined for $N_e<9$, and overdetermined for $N_e\ge 9$.
%Note that the results of Fig.~\ref{Figure4} validates our hypothesis
%regarding why $\epsilon_m$ is above $(\epsilon_m)_o$ for
%$N_e=5$ to $N_e=6$, and why $\epsilon_k$ is above $(\epsilon_k)_o$ for
%$N_e=10$ to $N_e=14$.  
%As stated in Part I, the condition number is an important
%factor in estimating where numerical instabilities may arise when
%solving the least squares problem (see part I).  
We point out that
a large condition number implies that the {\em upper bound} on the numerical
error may be large, but does not necessarily imply that the error
itself is large.  On the other hand, a very small condition number
does imply that the numerical error itself will also be small.
The variation in the instability number as a function of $N_e$ shown
in Fig.~\ref{Figure5} does in fact track fairly well the 
accuracy of the computed solutions
shown in Fig.~\ref{Figure4}, in the sense of the upper bound.
For this set of system parameters, a matrix with the logarithm
of the condition number
close to $10^{10}$ or higher has the potential to give rise
to numerical errors.

\begin{figure} 
\centering\scalebox{0.8}[1.0]{\includegraphics{massstiff2NODES.ps}}
\caption{The mass and stiffness
error parameters, $\epsilon_m$ and $\epsilon_k$, as a function
of $N_e$, for the
system parameters of Tables 1 and 2.  Lumped
elements are added to masses 3 and 21.
The horizontal lines represent the mass and stiffness
error parameters of the analytical model, $(\epsilon_m)_o$
and $(\epsilon_k)_o$.}\label{Figure4}
\end{figure}

\begin{figure} 
\centering\scalebox{0.8}[1.0]{\includegraphics{conditioncomparenew.ps}}
\caption{The instability measures for the
mass and stiffness updates as a function
of $N_e$, for the
system parameters of Tables 1 and 2.  Lumped
elements are added to masses 3 and 21.}\label{Figure5}
\end{figure}
%\item{Figure 5:  The instability measures for the
%mass and stiffness updates as a function
%of $N_e$, for the
%system parameters of Tables 1 and 2.  Lumped
%elements are added to masses 3 and 21.}

Figure~\ref{Figure6} shows the variations of the eigenvalue and modal
matrix error parameters, $\epsilon_{\lambda}$ and $\epsilon_{X}$,
as a function of $N_e$, for the case of adding 2 lumped masses
at nodes 3 and 21.  
%It is interesting to note that $\epsilon_{\lambda}$
%follows the peak patterns for both
%the mass and stiffness errors.  Namely, 
%$\epsilon_{\lambda}$ becomes large in the same regions
%in which $\epsilon_{m}$ and $\epsilon_k$ have peaks.
%seems almost to be an additive combination of the mass and stiffness errors.
%On the other hand, the 
We note that both the eigenvalue and
modal matrix error parameters increase approximately where
the mass matrix itself suffers from error.  In this
set of experiments, stiffness errors were not sufficiently 
pronounced to allow themselves to be evidenced in these plots.
%Thus, a small modal matrix error, $\epsilon_{X}$, implies
%that the updated mass matrix is accurate, but it
%may not guarantee that the updated stiffness matrix
%is good.

\begin{figure} 
\centering\scalebox{0.8}[1.0]{\includegraphics{EigvalModal2NODES.ps}}
\caption{The eigenvalue and modal matrix
error parameters, $\epsilon_\lambda$ and $\epsilon_{X}$, as a function
of $N_e$, for the
system parameters of Tables 1 and 2.  Lumped
elements are added to masses 3 and 21.
The horizontal lines represent the mass and stiffness
error parameters of the analytical model, $(\epsilon_\lambda)_o$
and $(\epsilon_{X})_o$.}\label{Figure6}
\end{figure}

%If one is only interested in making use of the modal matrix, 
%this may provide some
%justification for those updating approaches which assume that the
%stiffness matrix is exact.  
%********PHIL, what do you think of that last sentence?**********

%Thus, the condition number of the linear
%system used to perform the update depends on the
%number, location and magnitude of the added lumped
%elements, which in turn affects the quality of the
%update.





For the set of system parameters of Tables 1 and 2, 
we have seen that as few as 2 added lumped 
masses can be used to execute
the update, provided the masses are strategically placed.  
However, even when the updating schemes are numerically stable for
2 added masses, the updated mass and stiffness parameters
are not as accurate as those obtained when 5 added masses 
are used to perform the update.  
Thus, there appears to be a trade-off between numerical accuracy
and experimental effort.  When more lumped masses are
attached, more work and longer down-time will result; however,
the resulting updated model will be more accurate.  Conversely,
when fewer lumped masses are added, less work will be required but
at the expense of less accuracy.
To validate the above conjecture, lumped masses (of magnitude 0.2 kg)
are added to every single node.  Figures \ref{Figure7}
and \ref{Figure8} show the resulting
variations of $\epsilon_{m}$, $\epsilon_{k}$, $\epsilon_{\lambda}$
and $\epsilon_{X}$ as a function of $N_e$.
Note that for this case of 25 added lumped masses, 
the error parameters $\epsilon_{m}$, $\epsilon_{k}$,
$\epsilon_{\lambda}$ and $\epsilon_{X}$
decrease at a faster rate compared to those
of Figs.~\ref{Figure2} and \ref{Figure3}, 
implying that the updated system parameters
converge to the actual values more rapidly.  Thus, the above trade-off
hypothesis is validated.

\begin{figure} 
\centering\scalebox{0.8}[1.0]{\includegraphics{massstiffALLNODES.ps}}
\caption{The mass and stiffness
error parameters, $\epsilon_m$ and $\epsilon_k$, as a function
of $N_e$, for the
system parameters of Tables 1 and 2.  Lumped
elements are added to every mass of the system.
The horizontal lines represent the mass and stiffness
error parameters of the analytical model, $(\epsilon_m)_o$
and $(\epsilon_k)_o$.}\label{Figure7}
\end{figure}

\begin{figure} 
\centering\scalebox{0.8}[1.0]{\includegraphics{EigvalModalALLNODES.ps}}
\caption{The eigenvalue and modal matrix
error parameters, $\epsilon_{\lambda}$ and $\epsilon_{X}$, as a function
of $N_e$, for the
system parameters of Tables 1 and 2.  Lumped
elements are added to every mass of the system.
The horizontal lines represent the mass and stiffness
error parameters of the analytical model, $(\epsilon_\lambda)_o$
and $(\epsilon_X)_o$.}\label{Figure8}
\end{figure}

The results of Figs.~\ref{Figure2}, \ref{Figure4} and \ref{Figure7} indicate that
for the system of Fig.~\ref{Figure1},
fewer modes are required to update the mass
matrix than the stiffness matrix.  This can be explained by
noting that to render the least squares problem overdetermined
for the mass matrix, at least $\sqrt{N}$ measured modes are required, 
while to render the problem overdetermined for the stiffness
matrix, at least $\sqrt{3N-2}$ measured modes are needed.
%
% (see section 4.3).
%Finally, the results we presented illustrate the utility of the
%proposed mass and stiffness updating schemes.  The newly developed
%updating algorithms show considerable promise, and the authors
%plan to apply the proposed algorithms to update an analytical
%model using real data.




\section{Conclusion}

New mass and stiffness
updating algorithms are developed.  Using the
original test data and the newly acquired mass-modified
modes of vibration, the
mass and stiffness matrices of the analytical model
can be accurately corrected.
By manipulating the matrix equations in such a way
that the unknown correction mass and stiffness
matrices appear as column vectors, 
the connectivity information can be easily implemented, 
thus preserving the physical 
configuration of the system
and reducing the amount of computational effort required
to correct the analytical model.
In addition, the structure of the least
squares problems reveals
the minimum number of modes that would be needed
to measure in order to ensure a sufficiently accurate updated model.







When the set of measured modes is
complete, the updating routines return updated mass
and stiffness matrices that are nearly identical to those
of the actual system.  While the level of accuracy diminishes
as the number of measured modes is decreased, the resulting
updated model that is obtained is still more accurate than the
initial analytical model, even for a limited number of measured modes.  
To ensure the numerical stability
of the updating schemes, any masses
can be added to the actual system as long as the 
measured eigenvalues of the actual 
system and the mass-modified
systems are distinct.
This implies that the magnitude of the added masses can be
small, and the number of added masses can be
few.  The proposed mass and stiffness updating algorithms 
do require additional work and cause down-time,  
because the modes of vibration for 
the mass-modified system  
need to be measured.  The additional effort, however,
is a relatively small price to pay for the ability to
correct the analytical model accurately.

%\section{Acknowledgements}

\acks
We would like to thank the reviewers who gave a thorough and
careful reading to the original manuscript.  We would also
like to acknowledge the reviewer who suggested an alternative
approach to updating the stiffness matrix, which we have
incorporated into the revised paper.   Their comments 
are greatly appreciated and have helped to improve the quality of 
the paper.  L.G. de Pillis
thanks Argonne National Laboratory and Harvey Mudd College for
supporting this research.  Argonne support is under U.S. Department of
Energy Contract W-31-109-ENG-38.


\begin{thebibliography}{mmm.m}  
%


%\bibitem {Gysin}
%H. P. Gysin 1990
%{\it 8th IMAC}, 195--204.
%Comparison of Expansion Methods of FE Modelling Error 
%Localization.
%
%\bibitem {Imregun} M. Imregun and W. J. Visser 1991
%{\it Shock and Vibration Digest} {\bf 23}(1), 9--20.
%A Review of Model Updating Techniques.
%
%\bibitem {Zak}M. Zak 1983
%{\it ASCE Journal of Engineering Mechanics} {\bf 109}(2), 111--132. 
%Discrete Model Improvement by Eigenvector Updating.
%
\bibitem {Mottershead}
J.~E.~Mottershead and M.~I.~Friswell,
Model updating in structural dynamics:  A survey,
\textit{Journal of Sound and Vibration} 1993; \textbf{167}(2):347--375.
%
%\bibitem {Friswell}
%M. I. Friswell and J. E. Mottershead 1995
%{\it Finite Element Model Updating in Structural Dynamics}.
%Kluwer Academic Publishers, Norwell, Massachusetts.
%

\bibitem {Berman}
A.~Berman,
Mass matrix correction using an incomplete set of measured modes,
\textit{AIAA Journal} 1979; \textbf{17}(10):1147--1148.

\bibitem {Baruch}
M.~Baruch and I.~Y.~Bar Itzhack,
Optimal weighted orthogonalization of measured modes,
\textit{AIAA Journal} 1978; \textbf{16}(4):346--351.


\bibitem {Wei}
F.~S.~Wei,
Stiffness matrix correction from incomplete test data,
\textit{ AIAA Journal} 1980; \textbf{18}(10):1274--1275.

%
%\bibitem {BermanNagy}
%A. Berman and E. J. Nagy 1982
%{\it AIAA Journal} {\bf 21}(8), 1168--1173.
%Improvement of a Large Analytical Model Using Test Data.
%


\bibitem {Kabe} 
A.~M.~Kabe,
Stiffness matrix adjustment using mode data,
\textit{AIAA Journal} 1985; \textbf{23}(9):1431--1436.


\bibitem {Berman2}
A.~Berman,
Comment on optimal weighted orthogonalization of measured modes,
\textit{AIAA Journal} 1979; \textbf{17}(8):927--928.


\bibitem {ChenKuo}
J.~C.~Chen, C.~P.~Kuo, and J.~A.~Garba,
Direct structural parameter identification by modal test results,
\textit{24th AIAA/ASME/ASCE/AHS Structural Dynamics and Materials
Conference} 1983; 44--49.


%\bibitem {Nalitolela} 
%N.~Nalitolela, J.~E.~T.~Penny and M.~I.~Friswell,
%`Updating structural parameters of a finite element model
%by adding mass or stiffness to the system,'
%{\it 8th IMAC}, 836--842 (1990).


\bibitem {Gysin}
H.~P.~Gysin,
Comparison of expansion methods of fe modelling error localization, 
\textit{8th IMAC} 1990; 195--204.



\bibitem {Imregun} 
M.~Imregun and W.~J.~Visser,
A review of model updating techniques,
\textit{Shock and Vibration Digest}, 1991;
\textbf{23}(1):9--20.


\bibitem {Bathe}
K.~J.~Bathe, 
\textit{Finite Element Procedures in Engineering Analysis}.
Prentice-Hall, Inc.: New Jersey, 1982. 


\bibitem {Datta}
B.~N.~Datta,
\textit{Numerical Linear Algebra and Applications}.
Brooks/Cole Publishing Company: California, 1995.


\bibitem {GvL}
G.~H.~Golub and C.~F.~Van Loan,
\textit{Matrix Computations}, Second Edition.  
Johns Hopkins University Press: Baltimore, 1989.


%\bibitem {CourantHilbert}
%R. Courant and D. Hilbert 1953~
%{\it Methods of Mathematical Physics, Vol. I}.  
%Interscience Publishers, Inc., New York.



%\bibitem {Targoff} 
%W. P. Targoff 1976
%{\it AIAA Journal} {\bf 14}(2), 164--167.
%Orthogonality Check and Correction of Measured Modes.
%
%\bibitem {Wada} 
%B. K. Wada, C. P. Kuo and R. J. Glaser 1986
%{\it AIAA paper 86-0905}.
%Multiple Boundary Condition Tests (MBCT) for Verification of
%Large Space Structure.
%

%\bibitem {Lallement} 
%G. Lallement and S. Cogan 1992
%{\it 10th IMAC}, 487--493.
%Reconciliation between Measured and Calculated Dynamic
%Behaviours:  Enlargement of the Knowledge Space. 





\bibitem{gams}
R.~F.~Boisvert, S.~E~Howe, D.~K.~Kahaner and J.~L.~Springmann,
\textit{Guide to Available Mathematical Software}.
U.S. Department of Commerce and National Institute
of Standards and Technology, NISTIR 90-4237, 1990.


\bibitem {courant}
R.~Courant and D.~Hilbert,
\textit{Methods of Mathematical Physics}, Volume I.
Interscience Publishers, Inc.: New York, 1953.



\bibitem {meiro}
L.~Meirovitch, 
\textit{Computational Methods in Structural Dynamics}.  
Sijthoff \& Noordhoff: The Netherlands, 1980.




%\bibitem {GvL}
%G. H. Golub and C. F. Van Loan 1996
%{\it Matrix Computations}.  
%Johns Hopkins University Press, Baltimore.


%\bibitem {CourantHilbert}
%R. Courant and D. Hilbert 1953~
%{\it Methods of Mathematical Physics, Vol. I}.  
%Interscience Publishers, Inc., New York.




%\bibitem {Bathe}
%K. J. Bathe 1982 
%{\it Finite Element Procedures in Engineering Analysis}.  
%Prentice-Hall, Inc., New Jersey.
%
%\bibitem {Targoff} 
%W. P. Targoff 1976
%{\it AIAA Journal} {\bf 14}(2), 164--167.
%Orthogonality Check and Correction of Measured Modes.
%
%\bibitem {Wada} 
%B. K. Wada, C. P. Kuo and R. J. Glaser 1986
%{\it AIAA paper 86-0905}.
%Multiple Boundary Condition Tests (MBCT) for Verification of
%Large Space Structure.


%\bibitem {Lallement} 
%G. Lallement and S. Cogan 1992
%{\it 10th IMAC}, 487--493.
%Reconciliation between Measured and Calculated Dynamic
%Behaviours:  Enlargement of the Knowledge Space. 

%\bibitem {Nalitolela} 
%N. Nalitolela, J. E. T. Penny and M. I. Friswell 1990
%{\it 8th IMAC}, 836--842.
%Updating Structural Parameters of a Finite Element Model
%by Adding Mass or Stiffness to the System.

\end{thebibliography}


%\centerline{Appendix A}

\appendix

Because the modal matrices $[X]$ and $[X_a]$
are always of full rank, we claim that
matrix $[G]$ will also be 
of full rank.  Consider the following triple product which
has the same form as the left hand sides of Eqs.~(\ref{massnew1}) 
and (\ref{stiffnew2}): 
\begin{equation} 
[X]^T [\delta Y] [W] \label{general}
\end{equation}
The known matrices $[X]$ and $[W]$ are of dimension $N \times N_e$, 
and the unknown matrix $[\delta Y]$ is of dimension $N \times N$.
Let the $(i,j)$th element of $[X]$, $[W]$ and $[\delta Y]$ 
be represented by $x_{ij}$, $w_{ij}$ and $\delta y_{ij}$, respectively.  Let 
the $j$th column vector of $[X]$ and
$[W]$ be denoted by 
${\bf x}_j$ and ${\bf w}_j$ ($j=1,\cdots,N_e$), respectively.
%where both ${\bf x}_j$ and ${\bf w}_j$ are of length $N$.
Expanding and manipulating Eq.~(\ref{general}) so that the elements of the
unknown matrix appear as an unknown vector, we get
\begin{equation} \label{a53}
[G]\;\delta{\bf y}
\end{equation}
where $\delta{\bf y}$, of length $N^2$, is given by
\begin{equation} 
\delta{\bf y}=\left[
\delta y_{11} \cdots \delta y_{1N} \;|\;
\delta y_{21} \cdots \delta y_{2N} \;|\; \cdots \;|\;
\delta y_{N1} \cdots
\delta y_{NN} \right]^T
\end{equation}
Element by element calculations reveal that 
matrix $[G]$ is of the form
\begin{equation} \label{g}
[G] = \left[
\begin{tabular}{c}
$[G_1]$ \\
\vdots\\
$[G_i]$ \\
\vdots\\
$[G_{N_e}]$\\
\end{tabular} \right]
\end{equation}
where each block $[G_i]$, of dimension $N_e\times N^2$, is given by 
\begin{equation} \label{gi}
[G_i]=\left[
\begin{tabular}{ccccc}
$x_{1i}{\bf w}_1^T$ &$ x_{2i}{\bf w}_1^T $ & $ x_{3i}{\bf w}_1^T $ 
& $ \cdots $ & $ x_{Ni}{\bf w}_1^T $\\
$x_{1i}{\bf w}_2^T$ &$ x_{2i}{\bf w}_2^T $ & $ x_{3i}{\bf w}_2^T $ 
& $ \cdots $ & $ x_{Ni}{\bf w}_2^T $\\
\vdots  &\vdots &\vdots & &\vdots \\
$x_{1i}{\bf w}_{N_e}^T$ &$ x_{2i}{\bf w}_{N_e}^T $ & $ x_{3i}{\bf w}_{N_e}^T $ 
& $ \cdots $ & $ x_{Ni}{\bf w}_{N_e}^T $\\ 
\end{tabular} \right]
\end{equation}
We claim that $[G]$ has full rank, and we formalize this in 
the following theorem.

\noindent{\bf Theorem A:} 
{\it Given $N \times N_e$ real matrices $[X]$ and $[W]$ of full rank, if 
matrix $[G]$, of dimension $N_e^2 \times N^2$, is 
constructed from $[X]$ and $[W]$
as in Eqs.~(\ref{g}) and (\ref{gi}), then $[G]$ also has full rank.}

\noindent{\bf Proof:}
%
%Viewing the the first $N_e$ rows of $[G]$ as a set of vectors, we see that
%this set of vectors is necessarily linearly independent because of the given linear
%independence of the set of vectors ${\bf w}_j$ $(j=1,2,\cdots,N_e)$.
%
Because $[W]$ has full rank, the set of vectors ${\bf w}_j$ $(j=1,2,\cdots,N_e)$ is
linearly independent, and therefore block $[G_k]$ (which contains
rows $(k-1)N_e + 1$ through $kN_e$) is also linearly 
independent, for $k=1,\ldots,N_e$.
%That is, within each block there will be no linear dependence.  
We must next examine the possibility of linear dependence 
among rows in different blocks, and in particular, among rows 
numbered $k$, $k+N_e$, $k+2N_e$, $\ldots$, $k+(N_e-1)N_e$,
for $k=1,\ldots,N_e$.  

Without loss of generality, assume that there
is linear dependence to be found for $k=1$, which would imply that there exist
scalars $\alpha_j$ (for $j=1,\ldots,N_e-1$) so that we can write row $1$ of
block $1$ as a linear combination of the first
rows from subsequent blocks as follows:
\begin{eqnarray}
[x_{11}{\bf w}_1^T  \;\; x_{21}{\bf w}_1^T \;\; x_{31}{\bf w}_1^T \;\;\ldots 
\;\;x_{N1}{\bf w}_1^T] 
&=&
\alpha_1[ x_{12}{\bf w}_1^T \;\;  x_{22}{\bf w}_1^T  \;\; x_{32}{\bf w}_1^T  \;\;  
\ldots  \;\;  x_{N2}{\bf w}_1^T] \nonumber\\ 
&+&
\alpha_2[ x_{13}{\bf w}_1^T \;\;  x_{23}{\bf w}_1^T  \;\; x_{33}{\bf w}_1^T  \;\;  
\ldots  \;\;  x_{N3}{\bf w}_1^T] \nonumber\\ 
&+&
\cdots \\
&+&
\alpha_{N_e-1}[ 
x_{1N_e}{\bf w}_1^T \;\;  x_{2N_e}{\bf w}_1^T \;\;  x_{3N_e}{\bf w}_1^T \;\;  
\ldots \;\; x_{NN_e}{\bf w}_1^T] \nonumber
\end{eqnarray}
This is equivalent to writing
\begin{eqnarray}
x_{11}{\bf w}_1^T &=& (\alpha_1 x_{12} + \alpha_2 x_{13} + \cdots + 
\alpha_{N_e-1} x_{1N_e})
                    {\bf w}_1^T \nonumber\\ 
x_{21}{\bf w}_1^T &=& (\alpha_1 x_{22} + \alpha_2 x_{23} + \cdots + 
\alpha_{N_e-1} x_{2N_e})
                    {\bf w}_1^T \\
\vdots \nonumber\\
x_{N1}{\bf w}_1^T &=& (\alpha_1 x_{N2} + \alpha_2 x_{N3} + \cdots + 
\alpha_{N_e-1} x_{NN_e})
                    {\bf w}_1^T \nonumber
\end{eqnarray}
This implies that
\begin{equation} \label{59}
{\bf x}_1 = \alpha_1 {\bf x}_2 + \alpha_2 {\bf x}_3 
+ \cdots + \alpha_{N_e-1} {\bf x}_{N_e}
\end{equation}
Equation~(\ref{59}) says that the first column of $[X]$ can %${\bf x}_1$ can
be expressed as a linear combination of the last $N_e-1$ columns of $[X]$. This 
is clearly a contradiction, because $[X]$ is 
given to have full rank.  Therefore, our initial supposition that linear dependence
can be found among the rows of $[G]$ is shown to be impossible, and it follows that
$[G]$ must have full rank.


\end{document}  


\centerline{APPENDIX A}

The following first-order eigenvalue perturbation analysis is
well known and certainly not new.  The subsequent
derivation is shown primarily for the convenience and
benefit of the readers.

Consider the following generalized eigenvalue problem
\begin{equation} \label{16}
[{\cal K}_o]\y_{oj}=\rho_{oj}[{\cal M}_o]\y_{oj}
\end{equation}
where $[{\cal M}_o]$ and $[{\cal K}_o]$ represent the real, symmetric 
unperturbed system mass and stiffness matrices of size 
$N\times N$, whose $j$th eigenvalue and eigenvector 
are denoted by $\rho_{oj}$ and $\y_{oj}$, respectively.  
Assume the $N$ eigenvalues are all distinct
and the eigenvectors are properly normalized 
such that the following orthogonality conditions are met:
\begin{equation}\label{17}
\left\{ \matrix{
\y_{oj}^T[{\cal M}_o]\y_{oi} = \delta_i^j \cr
\y_{oj}^T[{\cal K}_o]\y_{oi} = \rho_{oj}\delta_i^j
\cr}
\right.  
\end{equation}
where $\delta_i^j$ is the Kronecker delta.


The system is now slightly modified or perturbed such that its
eigensolutions are obtained by solving
\begin{equation} \label{15}
[{\cal K}]\y_j=\rho_j[{\cal M}]\y_j
\end{equation}
where matrices $[{\cal K}]$ and $[{\cal M}]$ remain symmetric.
The perturbed system matrices are given by 
\begin{equation} \label{18}
\left\{ \matrix{
[{\cal K}]=[{\cal K}_o]+[\delta{\cal K}] + \cdots \cr
[{\cal M}]=[{\cal M}_o]+[\delta{\cal M}] + \cdots
\cr}
\right.  
\end{equation}
where $[\delta{\cal K}]$ and $[\delta{\cal M}]$ correspond to 
the symmetric first-order perturbation matrices.   
Assuming the modifications made are small,
then the eignesolution of Eq.~(\ref{16}) 
can be used to derive approximate expressions for the eigensolution 
of Eq.~(\ref{15}).  To this end, the eigenvalues and 
eigenvectors of Eq.~(\ref{15}) are assumed to
differ from those of Eq.~(\ref{16}) by some small 
perturbations as follows:
\begin{equation} \label{4}
\left\{ \matrix{
\rho_j=\rho_{oj}+\delta\rho_j + \cdots \cr
\y_j=\y_{oj}+\delta\y_j + \cdots
\cr}
\right.  
\end{equation}
where $(\rho_{oj},\y_{oj})$ are the $j$th unperturbed 
eigensolution, and $(\delta\rho_j,\delta\y_j)$ are the 
$j$th first-order perturbations.  The unperturbed 
eigenvectors are normalized according to
Eq.~(\ref{17}).  Thus, the unperturbed eigenvectors, 
the $\y_{oj}$'s, form a complete orthonormal set 
(with respect to $[{\cal M}_o]$) in the $N$-dimensional 
space, and any vector in 
that $N$-dimensional space may be expressed as a 
linearly combination of the $\y_{oj}$'s.
Hence, the $j$th first-order eigenvector 
perturbation can be written as
\begin{equation} \label{6}
\delta\y_j=\sum_{r=1}^N \epsilon_{rj}\y_{or}
\end{equation}
where the $\epsilon_{rj}$'s are the small coefficients 
to be determined.  

Substituting 
Eqs.~(\ref{4}) and (\ref{18}) into Eq.~(\ref{15}) and keeping only 
the first-order terms, we obtain
\begin{equation} \label{19}
[{\cal K}_o]\delta\y_j + [\delta{\cal K}]\y_{oj}=\rho_{oj}
[{\cal M}_o]\delta\y_j +
\rho_{oj}[\delta{\cal M}]\y_{oj} + \delta\rho_j[{\cal M}_o]\y_{oj}
\end{equation}
Substituting Eq.~(\ref{6}) into Eq.~(\ref{19}) yields:
\begin{equation} \label{20}
[{\cal K}_o]\sum_{r=1}^N\epsilon_{rj}\y_{or} + [\delta{\cal K}]\y_{oj}=
\rho_{oj}[{\cal M}_o]\sum_{r=1}^N \epsilon_{rj}\y_{or}+
\rho_{oj}[\delta{\cal M}]\y_{oj} + \delta\rho_j [{\cal M}_o]\y_{oj}
\end{equation}
Premultiplying Eq.~(\ref{20}) by $\y_{oj}^T$
and recalling the orthogonality conditions of Eq.~(\ref{17}), 
we have
\begin{equation} \label{21a}
\delta\rho_j=\y_{oj}^T \left( [\delta{\cal K}]-\rho_{oj}
[\delta{\cal M}] \right) \y_{oj}
\end{equation}
Thus, to the first-order, the $j$th perturbed 
eigenvalue is given by
\begin{equation} \label{24}
\rho_j=\rho_{oj}+\y_{oj}^T \left( [\delta{\cal K}]-\rho_{oj}
[\delta{\cal M}] \right)\y_{oj}
\end{equation}


\end{document}



