aspects of the present invention are disclosed in the following description and related figures directed to specific embodiments of the invention . those skilled in the art will recognize that alternate embodiments may be devised without departing from the spirit or the scope of the claims . additionally , well - known elements of exemplary embodiments of the invention will not be described in detail or will be omitted so as not to obscure the relevant details of the invention . as used herein , the word “ exemplary ” means “ serving as an example , instance or illustration .” the embodiments described herein are not limiting , but rather are exemplary only . it should be understood that the described embodiments are not necessarily to be construed as preferred or advantageous over other embodiments . moreover , the terms “ embodiments of the invention ”, “ embodiments ” or “ invention ” do not require that all embodiments of the invention include the discussed feature , advantage or mode of operation . further , many embodiments are described in terms of sequences of actions to be performed by , for example , elements of a computing device . it will be recognized that various actions described herein can be performed by specific circuits ( e . g ., application specific integrated circuits ( asics )), by program instructions being executed by one or more processors , or by a combination of both . additionally , these sequence of actions described herein can be considered to be embodied entirely within any form of computer readable storage medium having stored therein a corresponding set of computer instructions that upon execution would cause an associated processor to perform the functionality described herein . thus , the various aspects of the invention may be embodied in a number of different forms , all of which have been contemplated to be within the scope of the claimed subject matter . in addition , for each of the embodiments described herein , the corresponding form of any such embodiments may be described herein as , for example , “ logic configured to ” perform the described action . generally referring to exemplary fig1 - 6 , methods and systems for empirical modeling of time - varying , parameter - varying and nonlinear difference equations may be described . the methods and systems can be implemented and utilized to provide for a variety of results and which may be implemented efficiently . as shown in exemplary fig1 , a flow chart of a methodology for empirical modeling of time - varying , parameter varying and nonlinear difference equations according to an exemplary embodiment may be shown . here , in 102 , a set of time - varying , parameter varying and , if desired , nonlinear state space difference equations may be utilized . in 104 , the equations may then be expanded with respect to a chosen set of basis functions , for example nonlinear input - output equations may be expanded in polynomials in x t and u t . then , in 106 , difference equations may then be expressed as a linear time - invariant , for example , in terms of outputs y t and augmented inputs u t , which can include inputs u t and basis functions , for example polynomials , in inputs u t scheduling functions p t and states x t . exemplary fig2 may show an exemplary flow chart where a linear , parameter varying system of difference equations may be utilized . in this embodiment in 202 a set of linear parameter varying state space equations , such as those shown below as equation 1 and equation 2 , may be used . x t + 1 = a 0 x t + b 0 u t +[ a l ρ s ( 1 )+ . . . + a s ρ t ( s )] x t +[ b 1 ρ t ( 1 )+ . . . + b s σ t ( s )] u ( 1 ) y t = c 0 x t + d 0 u t +[ c 1 ρ t ( 1 )+ . . . + c s ρ ( s )] x t +[ d l ρ t ( 1 )+ . . . + d s ρ t ( s )] u t ( 2 ) then , in 204 , the state space difference equations may be expanded with respect to polynomials in the scheduling function p t , states x t and inputs u t , for example as shown in equations 3 and 4 below . x t + 1 = a 0 x t + b 0 u t +[ a 1 . . . a s ]( ρ t { circle around ( x )} x s )+[ b 1 . . . b s ]( ρ t { circle around ( x )} u t ) ( 3 ) y t = c 0 x t + d 0 u t +[ c 1 . . . c x ]( ρ t { circle around ( x )} x )+[ d 1 . . . d s ]( ρ t { circle around ( x )} u t ) ( 4 ) next , in 206 , the difference equations can be expressed in terms of original outputs y t and augmented inputs [ u t ,( p t { circle around ( x )} x t ), ( p t { circle around ( x )} u t )] that are functions of u t , x t and p t . difference equations can have linear time - invariant unknown ( a 0 , [ b 0 a ρ b ρ ], c 0 , [ d 0 c ρ d ρ ]) coefficients that can be estimated , and as shown in equations 5 and 6 below . then , as shown in 208 , the augmented inputs [ u t ,( ρ t { circle around ( x )} x t ), ( ρ t { circle around ( x )} u t )], and , in some exemplary embodiments , specifically ρ t { circle around ( x )} x t , can involve an unknown state x t vector , so iteration may be utilized or desired . thus , in such exemplary embodiments , iteration using an iterated algorithm , as described in further detail below , may be utilized . exemplary fig3 can show a flow chart of an iterated algorithm that may be implemented for iterated subspace identification . in this embodiment , and as seen in 302 , nonlinear difference equations can be expanded in additive basis functions and expressed in linear time - invariant form with augmented inputs u t . this can include , in some examples , nonlinear basis functions involving outputs , state and scheduling functions . then , in 304 , the state estimate { circumflex over ( x )} 6 [ 0 ] is unknown that is equivalent to ( a ρ c ρ ) term not in the lpv model , so the corresponding terms may be deleted from the set of augmented inputs ũ t . the iterated algorithm may then be implemented using augmented inputs as the inputs and can compute estimates { circumflex over ( θ )} [ k ] of the model parameters . then , using , for example , a kalman filter at the estimated parameter values { circumflex over ( θ )} [ 1 ] , the state estimates { circumflex over ( x )} t [ 1 ] can be computed along with the one - step prediction innovations . then the likelihood function can be evaluated . next , in 306 , an iteration k for k ≧ 2 may be made . here the state estimate { circumflex over ( x )} t [ k − 1 ] may be initialized for all t . the iterative algorithm may then be implemented using the augmented inputs as the inputs and can compute estimates { circumflex over ( θ )} [ k ] of the parameters . again , for example , using a kalman filter at the estimated parameter values { circumflex over ( θ )} [ k ] , the state estimates { circumflex over ( x )} t [ k ] may be computed . the one - step prediction innovations may also be made and the likelihood ρ ( y 1 : n | ũ 1 : n ;{ circumflex over ( θ )} [ k ] ) may also be evaluated . in 308 the convergence can be checked . here the change in the values of the log likelihood function and the state orders between iteration k − 1 and iteration k can be compared . if , in some examples , the state order is the same and the log likelihood function change is less than a chosen threshold ε . this threshold in many examples may be less than one , for example 0 . 01 , then the iterations may end or be stopped . otherwise , where the value is above a chosen threshold ε , step 306 above may be returned to and iteration k + 1 may be performed . following the performance of iteration k + 1 , the convergence may then be checked again . in another exemplary embodiment , a different approach may be taken to directly and simply obtain optimal or desired estimates of the unknown parameters for the case of autocorrelated errors and feedback in the system using , for example , subspace methods developed for linear time - invariant systems . this may be done by expressing the problem in a different form that can lead to a desire to iterate on the state estimate ; however , the number of iterations may be very low , and , to further simplify the system and its development , stochastic noise may be removed . for example , consider a linear system where the system matrices are time varying functions of a vector of scheduling parameters ρ t =( ρ t ( 1 ) ρ t ( 2 ) . . . ρ t ( s )) t of the forms as shown in the equations below : for affine dependence on the scheduling parameters , the state space matrices can have the form of the following equations 9 through 12 . in the above equations , it may be noted that where the matrices on the left hand side are expressed as linear combinations ( 1 ρ t ( 1 ) ρ t ( 2 ) . . . ρ t ( s )) t on the right hand side specified by the scheduling parameters ρ t , which may be called an affine form . in further discussion , the notation will be used for the system matrix a =[ a 0 a 1 . . . a s ]=[ a 0 a ρ ] and similarly for b , c , and d . in some further exemplary embodiments , system identification methods for the class of lpv systems can have a number of potential applications and economic value . such systems can include , but are not limited to , aerodynamic and fluid dynamic vehicles , for example aircraft and ships , automotive engine dynamics , turbine engine dynamics , chemical processes , for example stirred tank reactors and distillation columns , amongst others . one feature can be that at any given operating point pt the system dynamics can be described as a linear system . the scheduling parameters pt may be complex nonlinear functions of operating point variables , for example , but not limited to , speed , pressures , temperatures , fuel flows and the like , that may be known or accurately measured variables that characterize the system dynamics within possibly unknown constant matrices a , b , c and d . it may also be assumed that ρ t may be computable or determinable from the knowledge of any such operating point variables . for example lpv models of automotive engines can involve the lpv state space equations that explicitly express the elements of the vector ρ t as very complex nonlinear functions of various operating point variables . in some exemplary embodiments described herein , it may only be desired that the scheduling parameter ρ t may be available when the system identification computations are performed . this can be a relaxation of the real - time use or requirement for such applications as real - time control or filtering . to simplify the discussion , the lpv equations can be written in time - invariant form by associating the scheduling parameter ρ t with the inputs u t and states x t as x t + 1 = a 0 x t + b 0 u t +[ a 1 . . . a s ]( ρ t { circle around ( x )} x t )+[ b 1 . . . b s ]( ρ t { circle around ( x )} u t ) ( 13 ) y t = c 0 x t + d 0 u t +[ c 1 . . . c s ]( ρ t { circle around ( x )} x t )+[ d 1 . . . d s ]( ρ t { circle around ( x )} u t ) ( 14 ) here , { circle around ( x )} can denote the kronecker product m { circle around ( x )} n , defined for any matrices m and n as the partitioned matrix formed from blocks of i , j as ( m { circle around ( x )} n ) i , j = m i , j n with the i , j element of m denoted as m i , j . also , the notation [ m ; n ]=[ m t n t ] t can be used for stacking the vectors or matrices m and n . equations 13 and 14 above can then also be written as shown below in the formats of equations 15 and 16 . as discussed in more detail below , the above equations can be interpreted as a linear time - invariant ( lti ) system with nonlinear feedback of f t =[( ρ t { circle around ( x )} x t );( ρ t { circle around ( x )} u t )] where the states x t and inputs u t can be multiplied by the time varying scheduling parameters ρ t . the feedback f t inputs can now be considered as actual inputs to the lti system . as shown in further detail below , the matrices [ a ρ b ρ ; c ρ d ρ ] of the lpv system description can be the appropriate quantities to describe the lti feedback representation of the lpv system . further , the above equations may now be described as shown below in equations 17 and 18 . thus for measurements of outputs and inputs ũ t =[ u t t ( ρ t { circle around ( x )} x t ) t ( ρ t { circle around ( x )} u t ) t ] t the time - invariant matrices can be ( ã , { tilde over ( b )}, { tilde over ( c )}, { tilde over ( d )})=( a 0 , [ b 0 a ρ b ρ ], c 0 , [ b 0 c ρ d ρ ]) respectively . also , in situations where x t in ρ t { circle around ( x )} x t may not be a known or measured quantity , a prior estimate of x t may be available or utilized and iterations may be used to obtain a more accurate or desired estimate of x t . in still further exemplary embodiments , an lpv system can be expressed as a linear time - invariant system with nonlinear internal feedback that can involve the known parameter varying functions ρ t . ( see , section 2 . 1 of nonlinear system identification : a state - space approach , vincent verdult , 2002 , ph . d . thesis , university of twente , the netherlands , the contents of which are hereby incorporated by reference in their entirety ). in this exemplary embodiment the system matrices p i of rank r i may be factored for each i with 1 ≦ i ≦ s using a singular value decomposition , such as that shown in equation 19 . the quantities may then be defined as the following , shown in equations 20 through 23 . next , internal feedback in the lti system p 0 =[ a 0 b 0 ; c 0 d 0 ] may be considered with outputs z t = c z x t + d z u t and with nonlinear feedback from z t to f t = ρ t { circle around ( x )} z t entering the lti system p 0 state equations through the input matrices b f and d f . the state equations for x t + 1 and y t of the feedback system may be shown below in equations 24 through 27 . with respect to the above and referring now to exemplary fig4 , if it can be assumed that for any time t there is no effect of the feedback f t on the output z t as in the feedback structure shown in fig4 and in equation 26 , then there is no parameter dependence in the linear fractional transformation ( lft ) description ( see k . zhou , j . doyle , and k . glover ( 1996 ), robust and optimal control , prentice - hall , inc ., section 10 . 2 , the contents of which are hereby incorporated by reference in their entirety ). as shown in fig4 , there may be a system 400 having two boxes , box 404 which may be a memory - less nonlinear feedback system and box 402 that can be a linear time - invariant system . therefore , in some exemplary embodiments , the parameter dependence can become affine . exemplary fig4 may be a schematic diagram of equations 15 and 16 . the state equation 15 involves the upper boxes in 402 while the measurement equation 16 involves the lower boxes in 402 . δt 422 is a time delay of sample duration with the right hand side of equation 15 at 444 entering 422 and the left hand side equal to state x t + 1 at leaving . this is a recursion formula similar to equations 14 and 15 , so the time index can be changed from “ t ” to “ t + 1 ” for the figure before the start of the next iteration and continuing until entering boxes 420 , 430 and 410 . scheduling parameters ρ t 406 , inputs u t 408 and outputs y t 446 are variables the upper four boxes are multiplication from left to right by b 0 418 , a 0 420 , b ρ 414 , a ρ 416 , respectively . similarly , the lower boxes are multiplication with a replaced by c and b replaced by d , depicted as d ρ 424 , c ρ 426 , d 0 428 and c 0 430 . in the feedback box 404 , the kronecker products involving ρ t and successively x t and u t are formed in 410 and 412 respectively . the pairs of boxes aligned vertically are multiplying respectively from left to right the variables x t , u t , ρ t { circle around ( x )} x t and ρ t { circle around ( x )} u t . additionally , arrows shown as touching a wire can symbolize addition . further , δt 422 can represent a time delay block of duration δt that can act similar to a date line for this exemplary embodiment . therefore the arrows in exemplary fig4 indicate a time flow or an actual sequence of operations , the flow may start at 406 , 408 and the output of 422 and proceed through the diagram . upon reaching δt 422 , all operations for sample t have been performed and , upon crossing through δt 422 , sample time t + 1 may begin . thus , for example , upon leaving δt 422 , the same quantity is maintained , but all of the time labels can be changed to t + 1 throughout the process shown in exemplary fig4 . further , as shown below , this can be equivalent to the lpv form shown in equations 15 and 16 where the state equations can be linear in the scheduling parameter vector ρ . in further exemplary embodiments and now using a definition of feedback defined as f t = ρ t { circle around ( x )} z , as given by equations 26 and 27 , the state equations for x t + 1 and y t may be as shown in equations 28 and 29 below . then , if equation 19 is used to define the above factors , equations 28 and 29 may be the same as equations 15 and 16 . next , in a further exemplary embodiment , as the rank of p i may not be known , the outputs z t may be set as the states and inputs z t =[ x t t ; u t t ] t that may subsequently be fed back through the static nonlinearity so that f t =[ ρ t { circle around ( x )} x t ; ρ t { circle around ( x )} u t ]. then └ c z , i e z , i ┘=[ i dim x 0 ; 0i dim u ] and equations 30 and 31 may be set as below . thus , from the above , the lpv coefficient matrix p i =[ a i b i ; c i d i ] can be the regression matrix of the left hand side state equation variables ( x t + 1 ; y t ) on the vector of nonlinear terms [ ρ i , t x t ; ρ i , t u t ]. more generally , the lpv coefficient matrix p ρ =[ a ρ b ρ ; c ρ d ρ ]=[ b f , d f ] can be the regression matrix of the left hand side state equation variables ( x t + 1 ; y t ) on the vector of nonlinear feedback terms ( ρ t { circle around ( x )} x t ; ρ t { circle around ( x )} u t ). the lti nonlinear feedback representation can solve a major barrier to applying existing subspace identification algorithms to the identification of lpv systems and overcomes previous problems with exponentially growing numbers of nonlinear terms used in other methods . for example , the above lti nonlinear feedback representation can make it clear that nonlinear terms ( ρ t { circle around ( x )} x t ; ρ t { circle around ( x )} u t ) can be interpreted as inputs to an lti nonlinear feedback system . therefore it may be possible to directly estimate the matrices of the lti system state space equations using linear subspace methods that can be accurate for processes with inputs and feedback . this can directly involve the use of the outputs y t as well as augmented inputs [ u t ;( ρ t { circle around ( x )} x t ; ρ t { circle around ( x )} u t )] of the lti nonlinear feedback system . in another exemplary embodiment , lti system matrices and state vectors may be determined following the reduction of an lti subsystem of a nonlinear feedback system involving known scheduling functions and the state of the lti subsystem . this embodiment can involve taking the iterative determination of both the lti system state as well as the lti state space matrices describing the lti system . one example may be to consider the polynomial system as a linear system in x and u with additional additive input terms in the higher order product terms so the additional inputs are ρ t { circle around ( x )} x t and ρ t { circle around ( x )} u t . the scheduling variables p t are assumed to be available in real time as operating points or measured variables . if accurate estimates of the state x t are also available , then the problem could be only a direct application of the iterative algorithm for system identification . since the variables x t are not available until after the solution of the system identification , a different approach may be utilized . thus , in an exemplary first step , an initial estimate of the state vector may be made . here , system identification may be performed on the terms in the state equations involving the variables x t , u t and ρ t { circle around ( x )} u t but not the variables ρ t { circle around ( x )} x t . from this an approximation of the linear time invariant ( lti ) part of the system giving estimates of a 0 , b 0 , c 0 , d 0 , b ρ , and d ρ as well as estimates for the state vectors x 1 , n [ 1 ]=[ x t t x 2 t . . . x n t ] t may be obtained . then , in an exemplary second step , an iterate estimate of the state vectors may be made . here the state vector x 1 , n [ 1 ] can be used as an initial estimates for x t in the terms ρ t { circle around ( x )} x t in equations 15 and 16 . then the iterative algorithm can be applied to obtain an estimate of the system matrices a , b , c , and d and a refined estimate of x 1 , n [ 2 ] . further , this step may then be iterated until a desired convergence is achieved . exemplary steps one through three above may therefore work with only a few iterations . thus , in an exemplary manner in which the iterative algorithm can be used to address the previously known problem of lpv system identification . thus , the following is an exemplary discussion of using the iterative algorithm in directly identifying the coefficients f ij and h ij of the additive polynomial expansions of the nonlinear difference equation functions f ( x t , u t , v t ) and h ( x t , u t , v t ), respectively . this may be a very compact and parsimonious parameterization for such a nonlinear system . the iterative algorithm described herein for linear time - invariant systems can therefore be used with only a very modest increase in computational requirements . further , this exemplary use of the iterative algorithm to directly treat additive nonlinear terms in the state equations involving the state vector such as ρ t { circle around ( x )} x t as additional inputs . since the state x t may not be initially known , this may facilitate iteration to estimate the state sequence starting with only the linear time invariant ( lti ) part of the system , a 0 , b 0 , c 0 and d 0 . as seen in the above exemplary discussion following equations 13 and 14 , a linear parameter varying system that can be affine in the scheduling variables ρ t can be expressed in time invariant form involving the additional input variables ρ t { circle around ( x )} x t and ρ t { circle around ( x )} u t . note that this involves nonlinear functions of ρ t with x t and u t . the dynamic system can be linear time - invariant in these nonlinear functions of the variables . in a further exemplary embodiment , the effect of additional inputs can be traced through the iterative algorithm . the exemplary steps outlined below may further be reflected in the table shown in exemplary fig5 and the flow chart of exemplary fig6 . using elements 502 of exemplary fig5 and in step 602 of exemplary fig6 , an arx model may be fitted to the observations on iteration k with outputs y t and inputs [ u t t ( ρ t { circle around ( x )}{ circumflex over ( x )} t [ k − 1 ] ) t ( ρ t { circle around ( x )} u t ) t ] t , where for k = 1 the term ρ { circle around ( x )}{ circumflex over ( x )} t [ k − 1 ] may be absent from the inputs and further can be the equivalent to removing the term involving ( a ρ c ρ ) from the lpv model structure . the arx model fitting can have a linear regression problem that makes no prior assumptions on the arx order other than the maximum arx order considered . if the identified arx order is near a maximum considered , the maximum arx order considered can be doubled and the regression recomputed . in exemplary step 604 a corrected future can be computed . the effect of future inputs on future outputs can be determined using the arx model and subtracted from the outputs . the effect of this can be to compute the future outputs that could be obtained from the state at time t if there were no inputs in the future of time t . in exemplary step 606 , a canonical variate analysis ( cva ) can be made or computed . here , the cva between the past and the corrected future can be computed . again , the covariance matrices among and between the past and corrected future may also be computed . this may be similar to an svd on the joint past - future covariance matrix which is of the order of the covariance of the past to obtain the arx model . a result of this step is to obtain estimates of the states of the system called ‘ memory ’. in exemplary step 608 , a regression using the state equation may be performed . the ‘ memory ’ from step 606 can be used in the state equations as if it were data and resulting estimates of the state space matrices and covariance matrices of the noise processes can be obtained . these estimates can be asymptotically ml estimates of the state space model with no prior assumptions on the parameter values of the arx or ss model . thus , using the above - described methodology , the ml solution in the iterative algorithm can be obtain in one computation based on the assumed outputs and inputs in iteration k , as shown in 504 of exemplary fig5 , with no iteration on assumed parameter values . the iteration is the result of refinement of the state estimate in the nonlinear term ( ρ t { circle around ( x )}{ circumflex over ( x )} t [ k − 1 ] ) t that can be part of the assumed data in the iteration k ( 504 ). referring back to step 602 , the arx model fitting , the dimension of the augmented outputs is increased from dimu to dimua = dimu + dimp ( dimx + dimu ). to fit the arx model , the arx order lagp identified can be substantially higher due to the nonlinear input terms and depending on the statistically significant dynamics present among the output and augmented input variables . the computation can involve computation of an svd on the data covariance matrix that is dimension of the dimua * lagp where lagp is the maximum arx order considered . the computation that may be utilized for the svd is of the order of 60 *( dimua * lagp ) 3 , so the computation increases proportional to ( dimp ( dimx + dimu )/ dimu ) 3 . therefore , one consequence of augmenting the system inputs by the nonlinear terms ρ t { circle around ( x )} x t and ρ t { circle around ( x )} u t may be to increase the past by a factor dimp ( dimx + dimu )/ dimu , and to increase the computation by this factor cubed . depending on the particular dimensions of u t , x t , and ρ t , this can be very significant , however there is no exponential explosion in the number of terms or the computation time . the lpv subspace algorithm of this invention still corresponds to subspace system identification for a linear time - invariant system and , in addition , because of the nonlinearity of the terms [ u t t ( ρ t { circle around ( x )} x t ) t ( ρ t { circle around ( x )} u t ) t ] t involved the state estimates x t , iteration on the estimate of the system states until convergence can be desired or , in some alternatives , required . another factor to be considered is that the result of the cva in exemplary step 606 is the computation of an estimate , denoted as m t , of the state sequence x t . the symbol ‘ m ’ is used as in the word ‘ memory ’ to distinguish it from the actual state vector and various estimates of the state vector as used in the em algorithm discussed below . the estimate { circumflex over ( m )} t in combination with the maximization step of the em algorithm can be shown to be a maximum likelihood solution of the system identification problem for the case of a linear time - invariant system . in that case , the global ml solution can be obtained in only one iteration of the algorithm . this may be different with the em and gradient search methods that almost always utilize many iterations . the cva estimate m t may actually be an estimate of the state sequence for the system with parameters equal to the maximum likelihood estimates in the case of lti systems and large sample size . this is different from the usual concept of first obtaining the ml parameter estimates and then estimating the states using a kalman filter . not only is the optimal state sequence for the ml parameters obtained in a single iteration , the optimal state order may also be determined in the process . in the usual ml approach , it can be desired or , in some alternatives , required to obtain the ml estimates for each choice of state order and then proceed to hypothesis testing to determine the optimal state order . in another exemplary embodiment , the convergence of the iterative algorithm for the case of lpv may be described . in this embodiment , a substantially similar approach may be taken for other forms of nonlinear difference equations . for any iterative algorithm , two issues may typically arise : ( 1 ) does the algorithm always converge , and ( 2 ) at what rate does it converge . in the computational examples considered , the iteration k on the estimated states { circumflex over ( x )} t [ k ] can be very rapid , for example half a dozen steps . thus it may be shown herein that the iterative algorithm can be closely related to the class of em algorithms that can be shown to always converge , under an assumption on the lpv system stability . also , the rate of convergence can be computed to be geometric . this latter result is unique since the em algorithm typically makes rapid early progress but becomes quite slow in the end . the reason for the rapid terminal convergence of the lpv algorithm will be discussed in further detail . issues of initialization , stability and convergence will be elaborated below . additionally , the methods and systems described herein may be discussed in the context of the em algorithm as there can be some parallelism between the two . to show the convergence of the lpv algorithm , the development in gibson and ninness ( 2005 ), denoted as gn below and incorporated by reference herein in its entirety ( s . gibson and b . ninness , “ robust maximum - likelihood estimation of multivariable dynamic systems ,” automatica , vol . 41 , no . 5 , pp . 1667 - 1682 , 2005 ) will be discussed with various modifications made that are appropriate for the lpv algorithm . all equation numbers from gibson and ninness ( 2005 ) will include gn following the number in the paper gn , for example ( 23gn ). in this exemplary embodiment , the replacements that can be made in the gn discussion to obtain the lpv algorithm may be as follows : replace the lti state equations with the lpv state equations and , for the missing data , replace the state estimate from the kalman smoother with the ‘ memory ’ vector m t in the iterative algorithm . the consequence of this can be significant because for linear systems as in gn the iterative algorithm described herein can obtain the global ml parameter estimates in one step in large samples . on the other hand , for linear systems it may take the em algorithm many iterations to obtain the ml solution . replacing the lti model with the lpv model replaces the state space model equation 11gn by the lpv equations ( 15 ) and ( 16 ). this can produce a number of modifications in the equations of lemma 3 . 1 since u t and ( a , b , c , d ) in gn is replaced by [ u t t ( ρ t { circle around ( x )} x t ) t ( ρ t { circle around ( x )} u t ) t ] t and ( a 0 , [ b 0 a ρ b ρ ], c 0 , [ b 0 c ρ d ρ ], ) respectively . so the ‘ data ’ includes the vector p t { circle around ( x )} x t where the state vector x t is not available . to execute iteration k as shown in 504 of exemplary fig5 , a straight forward approach involving the em algorithm could be used . the ‘ missing data ’ would include state vector x t as the missing data . ( see s . gibson , a . wills , and b . ninness ( 2005 ), “ maximum likelihood parameter estimation of bilinear systems ”, ieee trans . automatic control , vol . 50 , no . 10 , pp . 1581 - 2005 , the contents of which are hereby incorporated by reference in their entirety ). then we could proceed for the case of a bilinear system that involves a term u t { circumflex over ( x )} x t that may be completely analogous to ρ t { circumflex over ( x )} x t for the lpv algorithm . but such an em approach can , on occasion , result in the typical slow convergence behavior of em algorithms near the maximum . therefore , instead the lpv algorithm may be used and can result in a rapid convergence . thus , in this exemplary embodiment , instead of specifying the ‘ missing data ’ as the estimate of the kalman smoother , the subspace approach can specify the cva state estimate m t or ‘ memory ’, as the missing data . the memory m t is the estimate of the state vector obtained by a canonical variate analysis between the corrected future and the past obtained in exemplary step 606 of the iterative algorithm using the input and output vectors specified in fig5 . this may be similar to a kalman filter state estimate at the global ml parameter estimates associated with the output and input data at iteration k rather than a kalman smoother state estimate at the last estimated parameter value . one difference can occur because the cva method expresses the likelihood function in terms of the corrected future conditioned on the past so that estimates of memory m t may depend only on the output and input data and their distribution depends on the ml estimates associated with the output and input data used in iteration k rather than smoothed estimates of the state . the actual conditional likelihood function p e ( ξ t / z t ) of ( 26gn ) with ξ t t =[ x t + 1 t y t t ] and z t t =[ x t t u t t ] is what can be involved in all of the em computations . this can be the same likelihood function involved in the exemplary step 608 of the cva algorithm of estimating the parameters of the state space equation as in lemma 3 . 3 of gn . a difference is that in the exemplary step 608 of the cva algorithm the expectation can be with respect to the true global ml estimates associated with the output and input data at iteration k whereas the gn estimate is an expectation with respect to the parameter value obtained in the previous iteration . further , the use of the lpv model and the choice of memory m t as the missing data in the expectation step of the algorithm can have the following consequence in gn . the basic theory in section 2 , the expectation - maximization ( em ) algorithm of gn , needs no modification . in section 3 of gn , the missing data is taken to be the cva state estimates { circumflex over ( m )} t [ k ] based on the input and output quantities for iteration k in 504 of fig5 . so the x t in equations ( 22gn ) through ( 28gn ) can be replaced by { circumflex over ( m )} t [ k ] . thus , lemma 3 . 1 of gn holds but also can achieve the global ml estimate associated with the input - output vectors of exemplary fig5 in one step . lemma 3 . 2 of gn can be replaced by the iterative algorithm to obtain the memory estimates { circumflex over ( m )} t [ k ] . lemma 3 . 3 of gn is the same result as obtained in the iterative algorithm . an additional step may be used to compute { circumflex over ( x )} t [ k ] from the estimates ⊖ [ k ] and the linear time - varying state equations given by the lpv state equations . this step may be desired to obtain the state estimate { circumflex over ( x )} t [ k ] for starting the next iteration k 1 . this step can further allow for the iterative algorithm to produce q ( θ , θ l )≧ q ( θ l , θ l ) with equality if and only if { circumflex over ( x )} t [ k ] ={ circumflex over ( x )} t [ k + 1 ] . then , in principle , the memory { circumflex over ( m )} t [ k ] could be used in place of the state estimate { circumflex over ( x )} t [ k ] . in fact , { circumflex over ( m )} t [ k ] projected on the recursive structure of the state equations in equations ( 43gn ) and ( 44gn ) can produce the ml state estimates asymptotically and the corresponding optimal filtered state estimates { circumflex over ( x )} t [ k ] . thus the use of { circumflex over ( m )} t [ k ] instead of { circumflex over ( x )} t [ k ] in the computation of ρ t { circle around ( x )}{ circumflex over ( x )} t [ k ] as part of the ‘ input data ’ can lead to essentially the same result except for some ‘ noise ’ in the results . but the computational noise can be avoided by the small additional computation of { circumflex over ( x )} t [ k ] . also , in some exemplary embodiments , there may not be a need for a ‘ robust ’ improvement to the iterative algorithm described herein since it has been developed using primarily singular value decomposition computations to be robust and demonstrated as such for more than a decade . an exception is possibly the computation of the filtered state estimate { circumflex over ( x )} t [ k ] that could be implemented using the square root methods of bierman ( g . j . bierman , factorization methods for discrete sequential estimation , academic press ( 1977 ); republished by dover , new york ( 2006 ), the contents of which are hereby incorporated by reference in their entirety ) if ill - conditioned lpv dynamic systems are to be solved to high precision . in still another exemplary embodiment , it may be demonstrated that the lpv system identification algorithm may converge at a geometric rate near the maximum likelihood solution . here the result for a linear system can be developed . the same approach may work for an lpv system , but the expressions below may be time dependent and can be of greater complexity . further , in order to simplify the derivation herein , the time invariant feedback from equations 17 and 18 can be considered with the substitution of notation ( ã ,{ tilde over ( b )},{ tilde over ( c )},{ tilde over ( d )}, ũ t , u t ) by ( a , b , c , d , u t , ũ t ). thus , equation 17 and equation 18 , with noise v t in innovation form may be written as equations 32 and 33 below . thus , from the above , it may be seen that for measurements of outputs y t and inputs u t =[ ũ t t ( ρ t { circle around ( x )} x t ) t ] t , the time - invariant matrices are ( a , b , c , d )=( a 0 ,{ b 0 a ρ b ρ ], c 0 ,[ b 0 c ρ d ρ ]), respectively . then , solving for v t in equation 33 and substituting in equation 32 produces equation 34 below . x t + 1 = ax t + bu t + k ( y t − cx t + du t )=( a − kc ) x t +( b − kd ) u t + ky t ( 34 ) next , through recursively substituting the right hand side of equation 34 for x t can provide equation 35 below . in equation 35 , j can contain the arx coefficients and p t can be the past output y t and inputs u t =[ ũ t ,( ρ t { circle around ( x )} x t ),( ρ t { circle around ( x )} ũ t )]. here ũ t can be the original inputs such that u t can include the nonlinear kronecker product terms . in a further exemplary embodiment , asymptotically for a large sample , j can be arbitrarily close to constant for a sufficiently large iteration k such that where l t can be time varying and can have the time varying scheduling parameters ρ t combined with terms l . also , if ρ t for all t is bounded , l t will be similarly bounded . now , when writing the state sequence in block vector form as x 1 : n i = vec [ x 1 : n i . . . x 1 i ] where the vec operation can stack the columns of a matrix starting with left hand columns on top , the above result can imply that where → can mean that l n − k extends to the right . using m to denote the upper triangular matrix can give the fundamental expression for the difference δx lag : n i between state sequences x lag : n i i at successive iterations k and k − 1 as equation 36 below . from equation 36 it may then be seen that the terminal convergence rate of the iteration in the lpv case can be governed by m , in particular the largest singular value of m . the various blocks of m can thus be computed by where the subscript ρ { circle around ( x )} x means to select the submatrix of b − kd with columns corresponding to the corresponding rows of ρ { circle around ( x )} x in ũ t . in some other exemplary embodiments , the convergence of the iterative linear subspace computation may be affected by the stability of the lpv difference equations and , more specifically , the stability of the lpv linear subspace system identification described herein . because a set of time - invariant linear state space difference equations may be stable if and only if all of the eigenvalues of the state transition matrix are stable , for example the eigenvalues are less than 1 . the lpv case is more complex , but for the purposes of this exemplary embodiment , the rate of growth or contraction per sample time can be given for each eigenvector component of the state vector x t by the respective eigenvalues of the lpv state transition matrix from equation 9 and now shown as equation 37 below . in equation 37 , the matrices a =( a 0 a ρ )=( a 0 a 1 . . . a s ) can be assumed to be unknown constant matrices . therefore , it is apparent that the transition matrix a ( ρ t ) may be a linear combination [ 1 ; ρ t ] of the matrices a i for 0 ≦ i ≦ s . therefore , for any choice of the matrices a i for 0 ≦ i ≦ s , there can be possible values of ρ t that could produce unstable eigenvalues at particular sample times t . if or when this occurs only sporadically and / or for only a limited number of consecutive observation times , no problems or errors may arise . in some other exemplary embodiments , for example in the k - th iterative computation of the subspace algorithm , only estimated values of the state sequence { circumflex over ( x )} t [ k ] and the matrices â i [ k ] for 0 ≦ i ≦ s may be used . therefore , if large errors in the estimates of these quantities are possible or acceptable , then there can be a greater potential for unstable behavior . for example , there may be areas of application of more significant importance , such as including identification of aircraft wing flutter dynamics where the flutter dynamics may be marginally stable or even unstable with the vibration being stabilized by active control feedback from sensors to wing control surfaces . then , in some other applications , it may be possible to guarantee that not any combination of scheduling parameters ρ t , values of matrices a i and uncertainty could produce instabilities . further to this , the stability of the transition matrix can provide a potential for issues and therefore further consideration could be desired . in some further exemplary embodiments , instabilities can produce periods where the predicted system response can rapidly grow very large . for example , when the eigenvalues are all bounded less than 1 , then the predicted system response can be bounded , whereas if estimation errors are large , for example , one of the eigenvalues is equal to 2 for a period of time , then the predicted response may double approximately every sample during that time . thus if 10 30 =( 10 3 ) 10 ≅( 2 10 ) 10 = 2 100 , then in as few as 100 samples , 30 digits of precision could be lost in the computation which could then provide meaningless results . therefore , per the above , there may be conditions under which there may be considerable loss of numerical accuracy that can be associated with periods where the lpv transition matrix is unstable , for example with extended intervals of time . further , the difficulty can lie in that the algorithm initialization and sample size since at the optimal solution with a large sample , the algorithm is stable and convergent . if it was possible to compute with infinite precision , then problems with illconditioning could be avoided ; however , with 15 or 30 decimal place accuracy , for example , some real data sets such as for the aircraft wing flutter , can benefit from further consideration . therefore , some exemplary embodiments may deal with manners of correcting for or otherwise lessening any undesired effects that may result from algorithm instability . for example , if the state sequence { circumflex over ( x )} t [ k ] is sufficiently close to the optimum as based upon the terminal convergence results described previously , the iterative algorithm may be stable provided it is assumed that the lpv system is stable . therefore , in some examples , large initial errors in the estimate { circumflex over ( x )} t [ 1 ] can lead to an unstable computation . in other examples , such as during time intervals when the scheduling parameter values cause significant instabilities in a ( ρ t ), some components of the f t | q t , the future outputs f t corrected for the effects of future inputs q t , may become quite large . therefore , in the computation of the covariance matrix of f t | q t , such large values can cause considerable illconditioning and loss in numerical precision . in examples where the iterative algorithm and some other subspace methods permit the arbitrary editing of the times where the instabilities occur , outlier editing of unstable regions may be performed . for example , in the arx covariance computation where the corrected future can be computed , the time intervals with significant instabilities can be determined and removed from the computation . in examples dealing with experimental design , if the trajectory of the operating point variables can be specified , then the scheduling parameters can be scheduled to enhance the system identification in several ways . an initial region that can avoid computational instabilities can be chosen to obtain sufficiently accurate initial estimates of the lpv parameters ( a , b , c , d ). this can then be used in other regions to initialize the state estimate of the algorithm with sufficient accuracy that computational instabilities will not be encountered . in the above examples , the removal of unstable outliers at each iteration can be the most general and robust procedure . as the estimated values of the state sequence { circumflex over ( x )} t [ k ] and lpv parameters improve with more iterations , the number of outliers can be expected to decrease until there is rapid terminal convergence . a counter example to this expectation is when the beginning of the scheduling parameters ρ t time history has little variation so that the lpv model for this part of the data is good for that portion of the data , but is a poor global model . then , in the later part of the time history , there can be consideration variation in ρ t such that unstable behavior may result . additionally , in many potential exemplary applications of the methods and systems described herein , it can be expected that the proposed algorithm can perform much better than existing methods that presently are not feasible on industrial problems . further , in many situations , it can be desired to design the experiment to obtain results of a desired fidelity for a specified global region of the operating space at as little cost in time and resources as possible . because the iterative algorithm &# 39 ; s linear subspace method is a maximum likelihood based procedure , designs can be developed for lti system identification . also , as it identifies a stochastic model with estimated disturbance models , including confidence bands on quantities such as dynamic frequency response functions , the required sample size and system input excitation can be developed with little prior information on the disturbance processes . in some further exemplary embodiments , the lpv methods and systems described herein may be extended to nonlinear systems . for example , it may be shown that a number of complex and nonlinear systems can be expressed in an approximate lpv form that can be sufficient for application of the lpv subspace system identification methods described herein . in one example , a general nonlinear , time varying , parameter varying dynamic system can be described by a system of nonlinear state equations , such as those shown in equations 38 and 39 . x t + 1 = f ( x t u t , ρ t , v t ) ( 38 ) y t = h ( x t , u t , ρ t , v t ) ( 39 ) in equations 38 and 39 , x t can be the state vector , u t can be the input vector , y t can be the output vector and v t can be a white noise measurement vector . in some exemplary embodiments , to deal with ‘ parameter varying ’ systems , the ‘ scheduling ’ variables ρ t that can be time - varying parameters can describe the present operating point of the system . very general classes of functions f (·) and h (·) can be represented by additive borel functions that need not be continuous . in a simplified manner , the case of functions admitting taylor expansion as in rugh section 6 . 3 ( w . j . rugh , nonlinear system theory : the volterra / wiener approach . baltimore , md . : johns hopkins univ . press , 1981 , the contents of which are hereby incorporated by reference in their entirety ), where ρ t and v t may be absent can provide the following , as shown in equations 40 and 41 . in equations 39 and 40 , the notation x t ( i ) can be defined recursively as x t ( i ) = x t { circle around ( x )} x t ( i − 1 ) and similarly for u t ( j ) . thus , equations 40 and 41 may be polynomial expansions of the nonlinear functions f (·) and h (·). note that the nonlinear equations may involve nonlinear functions of relatively simple form such as the approximating polynomial equations that involve only sums of products that are readily computed for various purposes . however for empirical estimation of the coefficients in the presence of autocorrelated errors , the problem can become difficult for low dimensions of y , u , and x , even using subspace methods . for subspace methods , the matrix dimensions can grow exponentially with the dimension of the ‘ past ’ that can be used . this can occur in expanding equation 40 by repeated substitution into x t [ k ] on the right hand side of the state equation 40 with x t on the left hand side of equation 40 . thus the entire right hand side of equation 40 with t replaced by t − 1 can be raised to the power lagp , the order of the past typically selected as the estimated arx model order . for a relatively low order past , and low dimensions of x t and u t , the number of additive terms can increase exponentially . however , in a further exemplary embodiment and following rugh the equations 39 and 40 can be converted through carleman bilinearization to bilinear vector differential equations in the state variable as shown in equation 42 : x t { circle around ( x )} =[ x t ( 1 ) ; x t ( 2 ) ; . . . x t ( i ) ]; ( 42 ) and the input power and products variables as shown in equation 43 . u t { circle around ( x )} =[ u t ( 1 ) ; u t ( 2 ) ; . . . u t ( j ) ]; ( 43 ) in a further exemplary embodiment , equation 40 , which expresses the state - affine form , can then be rewritten as equations 44 and 45 , below . x t { circle around ( x )} = a ( x t { circle around ( x )} { circle around ( x )} u t { circle around ( x )} )+ bu t { circle around ( x )} ( 44 ) x t { circle around ( x )} = c ( x t { circle around ( x )} { circle around ( x )} u t { circle around ( x )} )+ du t { circle around ( x )} ( 45 ) next , it may be assumed that the coefficient matrices ( a , b , c , d ) may be linear functions of scheduling parameters ρ t denoted as ( a ( ρ t ), b ( ρ t ), c ( ρ t ) d ( ρ t )) and in state - affine form as in equations 9 through 12 . the scheduling parameters ρ t may be nonlinear functions of the operating point or other known or accurately measured variables . for example , since the inputs u t { circle around ( x )} are multiplicative and can be assumed to be known in real time or accurately measured , they can be absorbed into the scheduling parameters ρ t , thereby possibly decreasing their dimension . the bilinear equations can then become those shown in equations 46 and 47 below . x t { circle around ( x )} = a ( ρ t ) x t { circle around ( x )} + b ( ρ t ) u t { circle around ( x )} ( 46 ) x t { circle around ( x )} = c ( ρ t ) x t { circle around ( x )} + d ( ρ t ) u t { circle around ( x )} ( 47 ) equations 46 and 47 are in explicitly lpv form . therefore , for functions f ( x t , u u , ρ t , v t ) and h ( x t , u u , ρ t , v t ) that may be sufficiently smooth and for which there may exist a well defined power series expansion in a neighborhood of the variable , there can exist an lpv approximation of the process . thus , with the absorption of the inputs u t { circle around ( x )} into the scheduling parameter ρ t , the equations may be linear time varying and , in some instances , may have a higher dimension ρ t . this can show the dual roles of inputs and scheduling parameters and how they may be interchangeable in a variety of manners and , for example , their difference may be their rate of variation with time . the knowledge of scheduling can therefore often be derived from the underlying physics , chemistry or other fundamental information . lpv models can therefore be considered , on occasion , as graybox models insofar as they may be able to incorporate considerable global information about the behavior of the process that can incorporate into the model how the process dynamics can change with operating point . the foregoing description and accompanying figures illustrate the principles , preferred embodiments and modes of operation of the invention . however , the invention should not be construed as being limited to the particular embodiments discussed above . additional variations of the embodiments discussed above will be appreciated by those skilled in the art . therefore , the above - described embodiments should be regarded as illustrative rather than restrictive . accordingly , it should be appreciated that variations to those embodiments can be made by those skilled in the art without departing from the scope of the invention as defined by the following claims .