设为首页收藏本站|繁體中文 快速切换版块

 找回密码
 立即加入
搜索
查看: 8562|回复: 4

时间序列达人们来讨论一下matlab中predict函数

[复制链接]
  • TA的每日心情
    开心
    2016-2-29 16:07
  • 签到天数: 1 天

    连续签到: 1 天

    [LV.1]初来乍到

    累计签到:1 天
    连续签到:1 天
    发表于 2010-9-28 11:13:27 | 显示全部楼层 |阅读模式

    马上加入,结交更多好友,共享更多资料,让你轻松玩转电力研学社区!

    您需要 登录 才可以下载或查看,没有账号?立即加入

    ×
    打开程序也看不懂,600多行代码,请达人们帮助解读一下,先谢谢各位仁兄了~~
    "真诚赞赏,手留余香"
    还没有人打赏,支持一下
    楼主热帖
    帖文化:【文明发帖 和谐互动】 社区精神:【创新、交流、互助、共享】

    该用户从未签到

    尚未签到

    发表于 2010-9-28 12:40:43 | 显示全部楼层
    贴代码上来哇
    "真诚赞赏,手留余香"
    还没有人打赏,支持一下
    帖文化:【文明发帖 和谐互动】 社区精神:【创新、交流、互助、共享】

    该用户从未签到

    尚未签到

    发表于 2010-9-28 13:03:22 | 显示全部楼层
    你想自己编一个predict函数?
    "真诚赞赏,手留余香"
    还没有人打赏,支持一下
    帖文化:【文明发帖 和谐互动】 社区精神:【创新、交流、互助、共享】

    该用户从未签到

    尚未签到

    发表于 2010-9-28 17:53:25 | 显示全部楼层
    恩,贴代码上来大家看看
    "真诚赞赏,手留余香"
    还没有人打赏,支持一下
    帖文化:【文明发帖 和谐互动】 社区精神:【创新、交流、互助、共享】
  • TA的每日心情
    开心
    2016-2-29 16:07
  • 签到天数: 1 天

    连续签到: 1 天

    [LV.1]初来乍到

    累计签到:1 天
    连续签到:1 天
     楼主| 发表于 2010-10-1 20:42:36 | 显示全部楼层
    回复 4# climber $ A7 B- w8 Y5 m. i

    ; `% ^# @1 d+ f% C* `) b" @0 Z) u, |1 H, q. E! l4 P' q: J( l
        谢谢大家的关注~~
    & l  M$ P) f" P2 Ffunction [yhat,xi,thpred]=predict(varargin)
    3 O8 K0 K7 M9 K* t; p" p  d%PREDICT Computes the k-step ahead prediction.4 R! y$ D: z2 O! `  V0 N5 e7 i/ d
    %   YP = PREDICT(MODEL,DATA,K)
    1 [) H; {8 b/ w, r$ N' B%! N- T/ h6 K& K4 x
    %   DATA: The output - input data as an IDDATA object, for which the  \0 f/ B" T  R7 N. }
    %   prediction is computed.
    2 X! C4 d. p) |%
    0 B; t# U' R$ _7 G" r- J%   MODEL: The model as any IDMODEL object, IDPOLY, IDSS, IDARX or IDGREY.) S- l1 D/ h2 r+ a& A2 j
    %
      y/ c) `% G4 v% e6 c6 g%   K: The prediction horizon. Old outputs up to time t-K are used to: S- M. k8 }! K1 A8 q
    %       predict the output at time t. All relevant inputs are used.
    ; }' |4 x3 {) X& j%       K = Inf gives a pure simulation of the system.(Default K=1).
    " l( k- I8 i) S0 r+ a' j& j%   YP: The resulting predicted output as an IDDATA object. If DATA- I( o) H: P: d, a
    %       contains multiple experiments, so will YP.# N- O- B) L2 f' G2 g
    %3 S3 O2 ?# N8 K2 N+ F& G5 \2 J1 X4 S; j
    %   YP = PREDICT(MODEL,DATA,K,INIT)  or
    . N2 {' C& G* W1 G%   YP = PREDICT(MODEL,DATA,K,'InitialState',INIT) allows the choice of
    ; l+ x) A; ]8 A/ A" q%   initial state vector:% ^. m4 h9 G+ z% y/ M
    %   INIT: The initialization strategy: one of
    / f( Q6 p/ k  `$ m%      'e': Estimate initial state so that the norm of the' D$ F. [' d+ B3 S" @+ {1 [; {
    %           prediction errors is minimized.+ C6 U# ?; T) w8 o3 U. _" s& o
    %           This state is returned as X0e, see below. For multiexperiment% J8 p2 d6 D1 u" v, ^- w- e9 b. G
    %           DATA, X0e is a matrix whose columns contain the initial states  y4 M8 U: S0 i
    %           for each experiment.
    6 N9 i  i& b: l( ~0 y# d. [%      'z': Take the initial state as zero
    4 V* Z' |) X& n7 h/ ^7 b' ^%      'm': Use the model's internal initial state.
    % z' _4 ~$ [" h+ z%      'd': Same as 'e', but if Model.InputDelay is non-zero, these delays
    2 E* m% [( y* p1 ^1 q+ N! R%            are first converted to explicit model delays, so that the are2 Z2 I$ V1 Y: Y; z2 [
    %            contained in X0 for the calculation of YP.- p( D! O! @* G; Z& }! `
    %       X0: a column vector of appropriate length to be used as initial value.
    # I5 c% J0 d3 o) l9 s' x5 q%           For multiexperiment DATA, X0 may be a matrix whose columns give
    ! F: n4 k/ G1 W%           different initial states for each experiment.
    # D+ y) L& x' H# V+ X( l6 ~9 [%   If INIT is not specified, Model.InitialState is used, so that+ P) Z5 k7 \) x* h7 A: X4 V' U5 h
    %      'Estimate', 'Backcast' and 'Auto' gives an estimated initial state," K0 m1 X2 e. Q5 L
    %      while 'Zero' gives 'z' and 'Fixed' gives 'm'. If Model is an IDARX8 P0 B6 n2 T' u
    %      model the default initial state is 'z'.
    # f9 {$ ~: v1 z& Z" M: ?%
    % h& c3 D) z+ |  |) h%   With [YP,X0e,MPRED] = PREDICT(MODEL,DATA,K) the initial state(s) and the
    * U' k5 D" N. C%   predictor MPRED are returned.  Note that if MODEL is continuous time, X0e" |' h. ^: F+ I4 y
    %   is returned as the states of this model. These may differ from the initial  ~: k5 \4 j; n
    %   states of the discrete time model that has the sampling interval(s) of
    # L5 z6 O1 q! e- O%   DATA. Also when INIT = 'd' only the states of MODEL are returned in* m3 y5 p4 V8 U
    %   X0e. (To obtain the full X0, apply first INPD2NK to a discrete time7 j2 c5 j7 J- _7 G5 E
    %   version of the model.)
    * u* n; o/ x' ^+ o8 Q%   MPRED is a cell array of IDPOLY objects,& s9 h5 Z$ G4 r) M( w. c
    %   such that MPRED{ky} is the predictor for the ky:th output. The matching
    & I' b: V4 }- T' g- A%   the channels of MPRED with data follows from its InputNames.  O  l. M) M' v  u( @3 w
    %   See also COMPARE and IDMODEL/SIM.+ l8 ]" K' o8 C5 d8 ]0 ~
    %   L. Ljung 10-1-89,9-9-94
    ; `7 J2 o$ h! m' j%   Copyright 1986-2008 The MathWorks, Inc.3 a8 O0 N0 H* `$ n/ b
    %   $Revision: 1.1.8.9.2.1 $  $Date: 2008/08/01 15:32:39 $
    5 c6 T7 D& {% W# _
    . Z6 k9 G; v3 m( J2 |/ e% First find if there is any pair 'InitialState', init in the argument! G; p9 _/ A9 @1 n
    % list:
    5 r3 \. Y2 ?, Z$ y. Jnr = find(strncmpi(varargin,'in',2));6 A9 w  t+ K1 K! X
    init =[]; xi = [];) s* @5 i% j2 `6 x" {/ ^3 h9 V
    if ~isempty(nr)
    ' {; E& i. K. O5 Z" {. A    if length(varargin)<nr+14 C  B- m! K  W
            error('Ident:general:optionsValuePair',...1 Y4 d5 O2 S  U, \
                    'The "%s" option of the "%s" command must be specified using a name/value pair. Type "help %s" for more information.',...
    - W, o, X) g2 u, W$ J                'InitialState','predict','predict')
    5 l# z  u) v8 @5 z9 a    end
    8 G" }3 n* u# i/ {    init = varargin{nr+1};9 S" ]7 F: b6 u, ]
        if ~isa(init,'char') && ~isa(init,'double')
    2 F0 B% N" `8 W. Y% B        error('Ident:analysis:IniSpec2',...
    ; W8 B" E1 [6 Y1 [* ?) _            'The "InitialState" option of the "%s" command must be set to a valid character or a vector/matrix. Type "help %s" for more information.','predict','predict')4 n6 i; [- `$ z8 N  {# g
        end6 F: k3 l: v! y. t
        varargin(nr:nr+1)=[];%6 W+ i0 |0 Z# x: e. H- T8 t
    end" n; M8 @& o* I& J) Y" e' I
    if length(varargin)<2. I( ]) j; I+ k+ N! X0 ?$ ]$ X% t
        disp('Usage: YP = PREDICT(MODEL,DATA)'), o6 T8 r" S" s6 ?) M
        disp('       YP = PREDICT(MODEL,DATA,M)')% D8 A% [. a- X- N. p
        yhat = []; xi=[];* ]2 }: t( p' d/ W6 e7 Q2 G$ P, Y0 @
        return; O# j6 c" p; Y3 I% X
    end# D2 _5 }3 ?  x- x1 C
    data = varargin{1};1 y, c5 d  i0 ^! p, q" u3 \3 _
    theta = varargin{2};# N5 A5 Q& g2 K+ l9 D: s
    if length(varargin)>2
    ( t1 e0 g; @  Z# A/ t% R/ H    m = varargin{3};9 _) J' a9 y  c1 L8 x! c
        if ~isnumeric(m) || ~isscalar(m) || ~isreal(m) || round(m)~=m || m<1% K) L7 x7 L! h  @
            error('Ident:analysis:predictInvalidHorizon',..." k; I- u5 b" {4 Z3 y4 [7 R+ w
                'In the "predict(MODEL,DATA,K)" command, the prediction horizon K must be a positive integer.')" O# ~/ J/ Z% ^8 Q+ ?% O
        end   
    3 t* _: b+ |, lelse
    ; U3 p5 A+ h+ J4 w* L" n, {    m = 1;( L1 q# L+ V+ }
    end
    3 }& L6 b& Y. v0 _# V: vif length(varargin)>3 && isempty(init): T9 f, ~: X, o6 j' r% Y
        init = varargin{4};2 e0 T6 z; N, T/ ?
    end& |& u7 |5 C( C1 @
    thpred = [];
    : |) ]! p& B+ i- f, o! s% ~3 \: |yhat = [];
    ) ]0 Z# @' a/ z' C, ?( N5 tif isa(data,'idmodel') % Forgive order) r: W: G8 p1 @/ ^
        data1 = theta;! Q& u2 I; n5 Z
        theta = data;
    7 F3 y9 o$ z7 H+ G( i    data = data1;- v5 Y, [: H( t  ?# C' X$ t. K
    end
    $ y& b: V/ g) {/ U5 F3 ~; enxorig = length(ssdata(theta));4 c0 Y7 J4 k( H! o+ Q3 N
    [nym,num] = size(theta);  ^! r! C9 [, V$ k" D
    if isa(data,'frd') || isa(data,'idfrd')! b0 a7 V. M* M/ D8 G: I
        error('Ident:analysis:predictIDFRD',.../ G4 S* \2 R, ^+ U. k  |; Q  Z
            'The "predict" command is not applicable to frequency response data (IDFRD object).')
    ; \& @' x2 }) b% G/ H( H  J  Xend
    " N3 ~% f# h& tif isa(data,'iddata')
    " V) U/ K* q. ~2 M4 T6 h    [Ncap,ny,nu] = size(data);$ ~8 x5 Y, t7 D& Z! J+ b
        if ny~=nym || nu~=num1 [. F% Q& a% {4 O; P6 r
            error('Ident:general:modelDataDimMismatch',...
    0 Z1 V! y4 k3 t0 {            'The number of inputs and outputs of the model must match that of the data.')
    . m$ b, p* ?" V6 v( i% z5 r0 B    end) i% t" {+ \, N* _( M
    end
    - R. A* }1 d% B5 j- m9 Tif isempty(m), m=1; end/ e+ ?" I8 _' s2 X
    if ~isinf(m),% T- K" H; Z1 z' C
        if m<1 || m~=floor(m)5 ^; I" u9 C' y
            error('Ident:analysis:predictInvalidHorizon',...
    $ t5 i! _3 {& B* d            'In the "predict(MODEL,DATA,K)" command, the prediction horizon K must be a positive integer.')* T" |- @' o7 N  D. p5 e
        end
    : w" |* u6 T9 U( Gend
    + X  ~* V9 B( w) @! H5 m3 Q5 ]% First deal with the special case of continuous time data and continuous
    3 M$ z) s' q& B3 l8 H  w% model:' K7 _9 P# }6 Q
    if isa(data,'iddata') && any(cell2mat(pvget(data,'Ts'))==0)% e- w) z+ o& \
        if ~isinf(m)
    + V, J( s# g9 M1 I+ Z        error('Ident:analysis:predictCTData',...
    0 O5 x3 r+ s  J) r+ G" E) p1 N; v; I            'The "predict" command cannot be used with continuous time data.')
    1 v( P1 w6 A2 k    end2 w1 G, `- A5 j& P
        yhat = sim(data,theta);%没懂" j. G* T8 I! _, J1 o& r
        if nargout == 0
    $ Q6 k$ R, n7 V: g; U2 i7 O/ e        utidplot(theta,yhat,'Simulated')
    " K: }" X: W8 z" p        clear yhat1 |, [+ V+ i# R
        end) j$ U! F( |$ I' z
        return; ~0 d8 o: E/ y3 H
    end
    ! T- }$ s% x$ V# tif isempty(init)
    ( m( Q& I- I4 R( y    if isa(theta,'idarx')
    2 Z2 W; R- X: W        init='z';1 t- t( k4 `5 j) l) G  C
        else7 H( S& v% e* U
            init=pvget(theta,'InitialState');
    + _8 o8 x" ~) f1 T) o    end$ q, g7 H2 e% F8 }- P
        init=lower(init(1));
    # q- e, |6 H+ Y4 U& g    if init=='d'! m' B: k( {! o! q9 H) g$ X. O
            inpd= pvget(theta,'InputDelay');
    5 |4 o9 [" ~! k- n5 A. q2 ?        if norm(inpd)==00 r9 L+ ~$ D' ~$ J
                init = 'e';
    4 q; c: P+ h$ f. l7 n2 V* U2 q  C3 G        end( a2 Q) l+ ?! p: j7 C3 m0 n
        end4 B7 J( E9 o6 S# {2 ?. |3 k  ?
        if init=='a' || init=='b', z; y, b! o6 i. q! r3 w; `
            init = 'e';
    2 A* V8 C3 q" S, w6 c& \    elseif init=='f'
    ' ?% Z, ^$ ~5 Z" t" @$ M, J9 @        init = 'm';* t) Y" V8 x" M
        end9 G5 j4 }: r9 b% p1 Y/ U8 z! X
    elseif ischar(init)
    . n* p, v& t' w7 \( L6 S    init=lower(init(1));3 X, W8 o  C# F* C" L1 S
        if ~any(strcmp(init,{'m','z','e','d'}))
    6 I# W7 w8 K  |' P        error('Ident:analysis:predictIni3',...* R# W+ c) `& N2 ^3 G' L0 n! x  l/ O) h
                'The "InitialState" option of the "idmodel/predict" command must be set to ''E(stimate)'', ''Z(ero)'', ''M(odel)'', ''D(elayconvert)'', or a vector.')
    9 l- b( A% P+ ^9 J: R: Z7 Q    end
    : E, f# q5 k5 d& C. \6 j( w9 C9 ~    if init=='m' && isa(theta,'idpoly')
    + p1 I! ?6 D( L8 @2 d        warning('Ident:analysis:predictIni4',...
    - r/ S( b# \# A# t2 w0 {            'Setting the "InitialState" option to ''m'' in the "predict" command is same as setting it to ''zero'' for IDPOLY models.')
    $ W  @8 n  t4 ]1 A& l( j( M    end
    9 {* @+ m6 \1 r2 x+ @. Tend. x3 w, \+ z+ K5 s1 Y+ m
    %% Now init is either a vector or the values 'e', 'm', 'z' 'd'
    6 O! v$ Q: D! }) w4 Iif isa(data,'iddata');
    6 ?6 ~' r$ S$ J2 S    if strcmp(pvget(data,'Domain'),'Frequency')
    $ w3 a& J3 {( X2 O/ U2 U        if m<inf
    ! Z, t/ S+ I, s7 c8 L5 N4 m            error('Ident:analysis:predictFreqDataFiniteK',...4 P+ b! K  {# u. W  R! L& a5 y0 j
                    'Prediction with finite horizon is not applicable to frequency domain data.')) t7 ]' K4 X. w  Y; R% y& z
            end
    / i/ N. ]; ~& H+ T    end3 t8 }( F& X( @* a3 y
    end
    # H3 r: ~( h* d  Bif isa(data,'iddata')
    2 H) |1 h# ], N6 @) U$ R    [uni,Tsd,inters] = dunique(data);
    . W; l; K" P: r1 m5 Q- z+ ^: X    iddflag = 1;, I- ]  t  {5 M6 h) D: D
    else
    + k- C9 N( S  c: H. C+ m    iddflag = 0;
    , z8 U! i4 y* [# h    uni = 1;
    & a/ [# Q0 J( `: _( t6 O2 \    Tsd = 1;) c- s: J8 {: x, T( t& f4 C
        inters = 'z';5 ]- o. v+ G! c% u* j/ U1 e6 N2 Q
    end
    # s, W; Q8 y7 y/ A%% Do away with continuous time right away:- t6 o9 }+ K/ n
    Ts = pvget(theta,'Ts');% n% M( w  T3 l! f! I& U( L7 R
    if Ts == 0
    % A7 M/ l2 s* L; L. p    theta = pvset(theta,'CovarianceMatrix',[]); % Less calculations- f+ B) R" p. ]
        ms = idss(theta);
    ! c  @( X, J' D( R) O/ }, l% }, V    nc = size(ms,'nx');; Y+ k. g, p# g* J" d8 u
        if ischar(init)
    / @4 B, ^( z$ C6 a" L# p# K6 L        if init=='m';
    " G2 I& k8 A1 n4 J* R+ t            init=pvget(ms,'X0');
    3 h" g) k5 ~" ~/ i% i& A( b        end
    0 W8 E  n3 y: ^# }4 F, e        if init=='z'
    6 r4 N# d& D0 u8 C9 j+ e            init=zeros(nc,1);
    ' Q7 p, F) |' |  L        end
      A' |) \+ p$ [. S    end9 Z7 C; R& S/ S/ a& x% n: }3 m2 t
        if ~isa(data,'iddata')
    # f; `$ e( I% k% y! V1 R9 L1 p        error('Ident:analysis:predictDataFormat',...
    8 @# j3 G- G, D( w9 @            ['For a continuous time model, data must be specified using an IDDATA object.',...
    : ?# I& D) ?4 G( [8 k            'Type "help iddata" for information on how to create the IDDATA object.'])$ _; c0 L$ A% E- n1 j
        end
    7 V% w) y$ x2 D9 X$ Q# t    Nex = size(data,'Ne');5 Z8 {1 r) M0 J6 N
        u = pvget(data,'InputData');
    ) }. O" a6 z" T& I* W( e    u1 = [];
    0 e! I! ~8 p& ]1 I    for kexp = 1:length(u);
    ! \* K  M7 L% J2 S  e. T0 O; X% \        u1 = [u1,u{kexp}(1,:).'];
    6 v7 C" J8 d3 Q) C) T& s, H    end
    # {# a3 Y& Z; n4 P2 j    nu = size(data,'nu');* ~+ K& J- E3 G$ `! O5 Z! D
        if uni: ]  T3 Q! g3 T$ W
            nc = size(ms,'Nx');
    6 Y) ^- ^6 `& a& u& ?: B        [md,Gll] = c2d(ms,Tsd,inters);
    . Y4 P4 n3 ?. \        if init=='d', md = inpd2nk(md); init='e';end" s9 {5 \! `& d3 {
            % First work on init/ u  f: N  X9 \; E3 F3 [' J
            if isa(init,'double')7 G' k! }* d* H6 Z; D, F2 r7 J( b
                [nxi,mexp]=size(init);
    % m9 ~, I% M* j) u( ?            if nxi~=nc
    9 C0 v5 W" c+ i: x3 q0 H! I- y                error('Ident:analysis:IniRows',...
    5 s1 b0 g( M6 C                    'The value of the "InitialState" option should have %d row(s).',nc)6 a" Z" g4 ~4 J
                end" G2 Y+ C$ v* k) U8 I
                if mexp==1 && Nex>1/ @+ l, k' n5 T! t
                    init = init*ones(1,Nex);/ q7 T% ]. m' I/ S. A$ A
                    mexp = Nex;3 e0 }& b- Z9 s
                end) i- J) A+ l. b+ r8 c
                if mexp~=Nex
    / O2 j/ o3 `3 z# J$ l                error('Ident:analysis:IniSize',...9 ^4 a2 x2 Q* o- l( r
                        'The "InitialState" option value must have either 1 or Ne columns, where Ne = no. of data experiments.');
    ' q* F- J% T1 R8 ^8 L1 x8 d            end1 }/ B  [8 J* ]: Z- R! s1 G( o9 k
                if ~isempty(Gll)2 w7 j3 T: o! o( x4 I
                    init = Gll*[init;u1];. b- h; d5 D7 T0 @6 D! h2 o; S
                end5 V$ m, l2 `$ [. J1 A% b# V
            end8 v- l, h9 M5 k% N8 t& ~
            if init == 'e'5 O2 I5 y  ]/ ^0 r- k, R7 P8 |# F9 D
                inite = x0iniest(md,data);
    # `, Z9 z& w. V0 s& p/ x        else: h" _- x3 m8 \5 O: a9 W
                inite = init;0 v2 H  ~4 Q% D5 D( c
            end) P7 v6 ^4 j# M5 S$ Y
            yhat = predict(md,data,m,inite); %Check XIC!
    / r' Z1 s7 h- q+ h& T! K        if ~isempty(Gll)
    & F- X. j3 u' H4 }+ j/ @+ E            xi = inite(1:nc,:)-Gll(1:nc,nc+1:nc+nu)*u1;
    " H% W/ F- }5 I8 d        else3 R( O& T7 i: c( r
                xi = inite(1:nc,:);6 V* T. R# b- V5 z. q
            end# G2 i5 g5 e# q0 ?; C2 ~5 ~; n
        else % Different sampling intervals  |3 l8 K7 O# t8 ]  @! d5 q. I
            Tsd = pvget(data,'Ts');; `& p6 k: u9 c" {" S
            ints = pvget(data,'InterSample');" T( E6 L7 `* e2 z) _
            %xic=zeros(nc,Nex);
    ) q! ]1 `- p6 B2 v& Q6 M0 k! k        for kexp = 1:Nex/ D; P! N! u% a# v$ Z9 l6 C
                [md,Gll] = c2d(theta,Tsd{kexp},ints{1,kexp});* B* {1 C' D5 y1 w4 I4 b, B
                if init == 'd', md = inpd2nk(md); init = 'e';end" {' Z4 T% O% U( u, ^  r/ Y0 m
                if isa(init,'double')& j/ S5 X8 m" G1 [2 |4 @# J* Y
                    if ~isempty(Gll)# k0 F' ^8 G7 ]( B: f7 }4 X
                        if size(init,2)==16 t, _6 P' R4 w
                            initk = Gll*[init;u1(:,kexp)];
    & K7 X" g1 q5 k; l; Q' `4 N  V                    else1 f# `8 E1 s+ J5 U9 m
                            if size(init,2)~=Nex
    8 Y& `8 @( k: g: ^5 v* m6 p                            error('Ident:analysis:IniSize',...
    . n4 o5 p( w6 y6 b+ ~. K                                'The "InitialState" option value must have either 1 or Ne columns, where Ne = no. of data experiments.');
    : G: O. o2 {2 w6 {) C/ M* O                        end5 n$ G" F6 L: A) k& B' ?
                            initk = Gll*[init(:,kexp);u1(:,kexp)];
    1 B2 a6 {' ]3 K" e0 J. r; q0 I) ]  `, s                    end% r5 z! D# @& _" O
                    end
    $ n3 {  ^! ?4 R7 M2 R. ?: |$ R7 R% M            else! x$ g* P& i8 O) L+ K% A, K
                    initk = init;
    ! d* Z& L! k3 c# ?" _" P2 R            end
    8 C0 g3 u/ C2 n; ]2 h0 L/ Q            if init == 'e'
    ) N: a0 u9 t% G: Y% q                initk = x0iniest(md,getexp(data,kexp));
    . Y( _, V% |; `) y5 E& p            end; f8 h$ p8 ^) e  q* Y+ C( Y
                if ~isempty(Gll)' x* _: m4 l/ U4 ?9 X0 ^! ?2 \
                    xi(:,kexp) = initk(1:nc)-Gll(1:nc,nc+1:nc+nu)*u1(:,kexp);
    0 W+ R, m+ I: w* R+ j7 y            else, a" _0 F1 ?8 S7 G; j/ ]
                    xi(:,kexp) = initk(1:nc);1 B0 \( j' \3 g2 B5 ?9 G
                end0 m, H4 L" a+ r. z( y& P+ l
                yhatk = predict(md,getexp(data,kexp),m,initk);
    0 S* `; n3 P$ W3 R( m            if isempty(yhat)  E, w, d- c+ S: n
                    yhat = yhatk;
    ; y1 {* z( ^4 N) b4 L            else
    4 V# H: K) |. a7 Z3 ]" |                yhat=merge(yhat,yhatk);
    & K" O. [& ]. {            end) A, z3 _. G: ]$ ^
            end
    ' ?+ w* j1 f" m6 h: `4 X    end
    / K8 |+ W, T' B8 {* f) a$ w6 l    if nargout == 0
    ' w4 B7 w: W- @% l        utidplot(md,yhat,'Predicted')' g  q* v+ p( S4 C- C* n7 y
            clear yhat, ?" J; @+ a6 F2 |
        end
    8 A0 U0 z! A" p    return
    , N% u% ~5 A0 [, C6 \else %Ts>0, G5 O7 f  k0 e  G5 d# |; z6 Q, ^: o
        if init=='d'. ]* U/ Z% {5 s( a+ h" p
            if any(pvget(theta,'InputDelay')>0)2 i7 B1 W$ b0 m1 f$ s& r
                theta = inpd2nk(theta);  u( |! c+ m" p
            end* J6 B) p+ |( {/ c7 x# _2 C* V
            init = 'e';& j7 r3 S' K* w  M0 t: Z  _- E
        end) {6 y1 E5 X- v
    end. j3 H+ d1 m0 H! U) v6 {. v% X( W
    if isinf(m)
    6 F0 w2 X+ }% G0 Q' ^5 v/ W    if ischar(init) && init=='e'7 z& {- x1 j  z3 |
            X0 = x0iniest(theta,data);$ J) k* O6 ~% `5 p0 L' A! J& ^
            init = X0;  % Note the difference with COMPARE, which fits
    " G% A9 }! O# q        % X0 with K = 0;
    9 X! _) Y+ l' t; h    end# `5 \2 D1 R! m! A, Y8 q0 C
        if ~isa(data,'iddata')) }% K0 C. b0 ^0 U! S
            data = data(:,nym+1:end);* x) n4 n1 S9 X! p2 E3 Z, h
        end+ A* H2 j! J0 v2 R0 s  D7 s9 t
        [lw0,id0] = lastwarn; was = warning('off','Ident:analysis:unstableSim');
    2 @* n4 g+ J$ T* d    try
    / m* \& }2 _( V8 g/ i) V4 b0 X7 L4 v        yhat = sim(theta,data,init);
    6 b+ u; h) I5 C. N  H1 V+ D0 n0 j. N    catch E" }/ h" J8 i- e
            warning(was)8 n  u8 ~2 I, g$ Z* m9 ?/ z
            lastwarn(lw0,id0);
    7 n" T3 d: U5 s3 Z        throw(E)$ G, @8 Y" ]& g2 P& ~! P% r, }
        end
    * B+ k# |3 [, L0 u0 x    warning(was)
    , n9 ^; M$ R6 l) o- [) _3 H, w    lastwarn(lw0,id0);# h5 [/ U  B3 D9 T
       
    5 K( @& \! t( s    thpred = theta;
    : F( ]* |% i$ [' K, w    if nargout == 0" q1 j( }1 `# o) X
            utidplot(theta,yhat,'Predicted')9 |. M5 a% D; Y" r# e: `+ n: ^$ U/ ?
            clear yhat4 z8 c* s- |  R; M! a& T) T4 M4 k7 s
        elseif nargout>1+ n+ Y3 V1 B. B0 ?
            xi = init;
    2 \" Z& t; j- t    end
    9 U$ I/ M" A0 O- B( U    return6 o) T5 R5 A; Q4 o& K
    end
    . P  {3 @! G' t" Y4 Q%Inpd = pvget(theta,'InputDelay');; ~. j2 d$ R. s' x1 ~
    nu = size(theta,'nu');) N1 u7 {' g( I- n  E
    if ~uni
    + s& v- @4 \+ Z4 i& e    Tsd = pvget(data,'Ts'); Tsd = Tsd{1};
    + T  n" o  F, G- }% C1 j1 K+ f; l    warning('Ident:estimation:nonUniqueDataTs',...
    , D7 \! y" y2 H7 h        ['The data set contains experiments with different sampling intervals.\n',...2 a8 c: M: P) m5 m* }6 @9 K- M
            'Sampling interval from the first experiment (=%g) will be used.'],Tsd)
    . R' o5 a6 x7 L" A$ c* v& O1 Qend$ Q7 l' F( D4 [4 x4 ^  }  d
    if iddflag && any(abs(Ts-Tsd)>1e4*eps)
    , q- T& g  U* A  k9 q, i    warning('Ident:analysis:dataModelTsMismatch',...$ o* L+ k& w) d& o' u3 t
                'The data and model sampling intervals are different. The data sampling interval will be used.')7 u/ I3 p* l. P1 ^! i
    end
    " u; Q" C, I9 t( Q: l5 Q; yInpd = pvget(theta,'InputDelay');
    6 B. k1 N! _7 |# jif init=='d'
    0 a% q+ Q  F$ |    init = 'e';
    9 @! i% U* n& a9 e! X9 b/ {8 G    theta = inpd2nk(theta);
    / }  c. {. n4 C6 P6 ?, t    Inpd = zeros(size(Inpd));
    . p8 L1 i" X0 r  Rend
    7 c7 Y* k: E3 _2 m, ]if isa(data,'iddata')7 F* i/ E/ b+ G5 h! O- \$ o& V3 _
        if isnan(data)
    3 a% X0 A* Z2 i        error('Ident:general:idprepMissingData',...
    4 N" u1 @$ w4 |6 C9 Q- ~& S1 V            'Data contains NaNs which represent missing values. Use the "misdata" command to fill in the missing values before using.');
    ( s; k# e1 i- d) t. n    end
    9 i$ a- z$ Q' W' K& b+ p5 g    data = nkshift(data,Inpd,'append');
    4 L* N& _1 U" R2 B    theta = pvset(theta,'InputDelay',zeros(nu,1));
    ' ]2 b! m/ A6 _    [ze,Ne,ny,nu,Ts,Name,Ncaps,errflag] = idprep(data,0,'dummy');4 @& s* I! l5 [* c! a
        %if ~isempty(errflag.message), error(errflag), end5 e% ?$ n9 n; w0 e- A
        if ~isempty(Name), data.Name = Name; end. b- A( m+ Y4 k+ Z
        iddatflag = 1;
    2 L5 i! p% x1 L" e" a- Xelse3 F# k0 T. E/ w
        if norm(Inpd)>eps
    5 |5 O: L/ y$ A- N' y/ n* V        if iscell(data)  2 D9 R4 G2 z" c) H' s: T3 j! n
                data = data{1};
    9 F! e9 I& e9 p0 }, y( s        end9 w" }6 T6 T0 ~! Q/ d, h+ C
            [Ncap,nudum] = size(data);
    / P# s6 q# D1 ?' k5 U        nk1 = Inpd;: \. R9 `# n- b- E
            ny = nudum - length(nk1);
    . I: c0 K0 h) J        Ncc = min([Ncap,Ncap+min(nk1)]);
    ; r3 R8 k) B% I( y: f4 u        for ku = 1:length(nk1)' F7 [7 q% ]. Z; c. Y) Q5 X" J
                u1 = data(max([nk1(:);0])-nk1(ku)+1:Ncc-nk1(ku),ny+ku);
    . C; Y( H) j4 f& E# S            newsamp = Ncap-length(u1);
    % P0 W  t" y7 v7 S, V            if nk1(ku)>0# v5 J. `$ w! b0 c
                    u1= [zeros(newsamp,1);u1];, c, z4 I, r1 @! S$ p
                else
    ; p9 d2 M7 y: s/ q                u1 = [u1;zeros(newsamp,1)];
    ; j0 q- J. _0 p9 I* h            end" g5 H8 _. ]) c  l" C, i
                data(:,ny+ku) = u1;/ j/ U# d' x' p0 E- C: J# Z3 ~
            end
    5 x: V/ E2 S0 C4 n' Z: b# e% Z) s; T4 S    end8 `; l& `+ d# a9 g; G. z
        iddatflag= 0;8 Z( g& N% `2 e/ d  n0 g
        if ~iscell(data), ze ={data};else ze = data;end7 j4 n# T* i6 G/ K9 i# d+ \
        Ne = length(ze);
    % r+ o) m; m; h1 p! @# j    nz = size(ze{1},2);
    . J7 H' _% J" b$ u: w2 I; n    for kexp = 1:Ne0 g; _. j6 M5 c# o0 x6 @3 K, o" L) f
            Ncaps(kexp) = size(ze{kexp},1);* n" l6 k$ B4 u0 Y  }- i6 o
        end
    6 I; Y: L0 q  N& b+ ~% c, C7 T$ N    %Ts =[];
    8 I9 m, t- M0 H5 J- J3 p& T! Cend% l  K7 |* f' {! R- b* T
    if m>=min(Ncaps). p/ S; t7 K( P4 a  |  J+ V
        error('Ident:analysis:predictLargeK',...0 f; r4 k/ @" W9 R
            'Prediction horizon must be smaller than the number of data samples.')6 S+ C0 P! ^3 e) X5 h5 w2 \
    end
    7 K. l# U  ?7 x$ P9 \1 b, W' k[A,B,C,D,K,X0]=ssdata(theta);6 I0 s$ \0 G& C7 p) i3 O
    [nyt,nx]=size(C);nut=size(B,2);
    + u# ]. `0 O9 C( F7 i7 X9 ~if iddatflag
    - a7 `: n) }. _5 E& p    if ny~=nyt || nu~=nut0 \- o. c( R( Q' y6 D
            error('Ident:general:modelDataDimMismatch',...
    4 @  x, x+ L* n' G1 @0 s* v            'The number of inputs and outputs of the model must match that of the data.'): j: ^' S( g. M1 u7 ~! S6 J4 j
        end7 O& T0 [9 ?1 H+ c! n
    else
    4 m; `9 Q0 f3 l! V3 p1 a$ W) u    if nz~=nyt+nut
    9 s5 ]+ H7 @" G        error('Ident:general:modelDataDimMismatch',...
    6 o$ b/ t7 s5 |            'The number of inputs and outputs of the model must match that of the data.')0 Y7 X+ y0 O+ |! X; g
        end5 Q% J- A1 D+ f4 B" p$ {
        ny = nyt; nu = nut;
    ) e* S& Z5 `; }end
    3 Y* ?, d5 ~0 J, Y" u9 Qif strcmp(init,'e')
    : X$ L  F' Q+ R- K    X0 = x0iniest(theta,ze);1 ]4 I' _+ o! n( [% {7 u6 d3 p
    elseif strcmp(init,'z'); O- M  Z& h$ B) z, g
        X0 = zeros(size(X0));1 ]2 F. z" ^! a$ J1 R9 k# f
    end
    ! @8 K2 i7 b' n  R) Cif ~ischar(init)* t8 n( V, a7 l( H/ ?
        X0 = init;! f2 {% a; h& Z  u# c5 k3 T
    end
    % q' E. v' [9 }8 J9 Eif nargout>1
    7 `& ]- a% t. U/ S8 @8 U8 W    xi = X0;
    # G3 r& F& S9 ]# p6 Kend
    5 b1 r5 G* b- L5 ?[xnr,xnc]= size(X0);( G8 j& a; _1 m" [( t+ j4 b$ D- ^* m
    if xnc~=1 && xnc~=Ne7 C# z( i! l: e* T: b8 P
        error('Ident:analysis:IniSize',...
    / Q; B5 q6 W  p$ M5 y# f: p        'The "InitialState" option value must have either 1 or Ne columns, where Ne = no. of data experiments.');' C" o* h4 |) ?0 ]9 R! B3 Z' @" d2 ?
    end
    4 |# j. f0 w) ~if xnr~=nx
    6 k) O9 B# ?  S6 B. G7 r9 z% X    error('Ident:analysis:IniRows',...
    ( k* _" e. [. J+ P/ F        'The value of the "InitialState" option should have %d row(s).',nx)
    3 G* g2 ?8 l& n2 _2 L& lend
    - r$ n$ x$ S5 P% B9 @! {if xnc==1 && Ne>1. [. H0 S! @- [0 x+ f" [
        X0= X0*ones(1,Ne);
    0 q5 @, ?1 _% N1 s& l  t( aend& n7 l: u4 C/ ^
    for kexp = 1:Ne
    % F7 N  R+ }, V8 N5 ^    z = ze{kexp};9 b0 V' ^; s$ u! x: H% M
        u = z(:,1+ny:end);8 j9 X- U' x$ L$ r4 d, P2 \  V6 V
        Ncap = Ncaps(kexp);
    8 ^: S4 u4 l7 ?  |+ H    if m==inf,3 z0 Q2 ~/ ?( J8 ~
            yhat=sim(theta,u,X0(:,kexp));  |* ~, {6 r* e: V1 u, Z
        else, d" D5 h3 c. `0 F: F. N
            x=ltitr(A-K*C,[K B-K*D], z, X0(:,kexp));$ X6 A( p" v# m. B' o
            if m==1,
    , `) r4 `: `  J0 [" {            yhat=(C*x.').';: e3 W/ ~5 h% z4 d7 P  X0 }/ o' {
                if ~isempty(D),yhat=yhat + (D*u.').';end
    7 W: I. g% {% V# K) O. U" {# r" g        else
    8 R  Q- H  J2 q: K! s8 m7 N- U$ M            F=D;Mm=eye(length(A));  J8 f$ l# E" V' F- |4 P& q
                for km=1:m-1
    7 s% z9 T9 l! K9 D4 k. X9 H                F=[F C*Mm*B];5 d; C- |7 j- u: W. k
                    Mm=A*Mm;
    9 f, a0 {4 ~6 N6 B            end" Y( U% m  Z/ q" s
                yhat=zeros(Ncap,ny);%corr 911111# U+ S' ?% \; N, z+ c  v
                for ky=1:ny
    + l  R; O8 e% e' e" v                for ku=1:nu
    ( z/ [* p, w) T# j8 r" \( ~                    yhat(:,ky)=yhat(:,ky)+filter(F(ky,ku:nu:m*nu),1,u(:,ku));
    5 S1 \( G8 t: n7 x9 [) k/ ~2 I4 t                end
    , C( l; O2 U) D6 ]  V            end
    7 I1 S0 o! C4 K) ?6 {* l            if isempty(yhat),yhat=zeros(Ncap,ny);end+ z5 i8 E) X7 |- V" T4 r2 y
                yhat(m:Ncap,:)=yhat(m:Ncap,:)+(C*Mm*x(1:Ncap-m+1,:).').';9 U, V  e# M: O
                if nu>0
    ! Y4 p, u* G0 }/ D$ B) ?                x=ltitr(A,B,u(1:m,:),X0(:,kexp));
    9 b  e2 i$ h! ^' P! j! C                yhat(1:m,:)=(C*x.').';
    % ]: R, }" j6 O  ~$ U& s            end9 K* z8 j! m* W. l+ k
                if ~isempty(D),yhat(1:m,:)=yhat(1:m,:) + (D*u(1:m,:).').';end2 P* R. Z  W' _
            end
    ) l$ p0 {# B: l( t% M, M    end
      j" P) ^  S6 Q5 n1 x6 n    yhatc{kexp} = yhat;! S! W2 B* |8 O7 Z, i6 D" h( I
    end
    5 x. s# F& e  G& `2 lif iddatflag: j: n0 D/ Z$ N0 U. ]; w7 r
        yhat = data;
    8 Z! U& D9 G2 Z- R# H    yhat = pvset(yhat,'OutputData',yhatc,'InputData',[]);
    : O' {; _0 F5 m1 {4 \else: N0 q1 W0 T/ Y  ^
        yhat = yhatc;/ _' e; M3 r  e* I
    end+ D8 {! V9 o' x- Q3 p
    if nargout >2
    - b+ E- S5 r& N! [. a% h    thpred = polypred(theta,m);
    % i% T' c7 B+ s+ m% Y  P) r: O: |end& O0 \+ g' ]/ {! d6 b$ F
    if nargout > 1
    % f# H" j  w: d$ `    try# a- V+ L1 p% o- k7 {5 O
            xi = xi(1:nxorig,:);7 i# @9 n# `/ m6 s; d8 W
        end5 `6 K0 G3 ^1 o3 x/ L- }5 q6 Y
    end
      _& }. F/ p  ~. v6 S; Qif (nargout == 0)
    1 [1 ]$ o" N  _. W  F    % Plot y and yhat., g$ d8 U& s# X3 |+ J8 Q( h
        utidplot(theta,yhat,'Predicted');
    : j5 ^  f6 t( q& e; ]    clear yhat x0;0 q* N: R/ t9 P/ \% \/ `. T% f5 O7 ?
    end
    4 N3 m. ^+ W; ?( P+ k  L2 p%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    ; [% w0 J0 a& v3 ?9 Q& r%%%%% LOCAL FUNCTIONS
    , V" d. }; J' S0 q: L%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%9 u5 T- }& n) y' k9 ~) q
    function xi = x0iniest(m,data)$ U5 Z7 R9 W8 M
    if isa(data,'iddata') && strcmpi(pvget(data,'Domain'),'frequency')
    2 y3 e* [% R, h+ o# y; c    [dum,xi] = pe_f(m,data);
    ( R# m. m! G  Q# x  `) s5 \    return9 g; b$ Q6 B4 Y2 j# W8 D3 Y) R
    end' P$ i# v2 e/ q) ?) T
    Inpd = pvget(m,'InputDelay');, T9 O& K9 P6 P$ Y9 j
    if isa(data,'iddata')6 r7 C- N; ]1 S' _1 r* D
        data = nkshift(data,Inpd,'append');
    , {  }; P) ]% }0 ^/ ]( S    [ze,Ne,ny,nu,Ts,Name,Ncaps,errflag] = idprep(data,0,'dummy');4 s8 {) O% K. N2 `3 M9 e
        %if ~isempty(errflag.message), error(errflag), end
    1 ?2 W3 d* L- Z  t: i    if ~isempty(Name)
    : H, b" G. C+ X# [  S3 q( S; g        data.Name = Name;3 n- ?, f. [$ Y* b  e% L$ B
        end2 M  r+ W/ {0 v/ l
    else
    ; Y& o$ x* [5 x( J9 N  D7 a) a    if iscell(data)
    0 E+ H, r! n; j: j# z! g- B( Q        ze = data;
    5 m2 P" [+ m- H! H    else
    : A+ O7 y5 w" `* o( [& U+ J( k2 D9 l        ze = {data};" d; z: V3 M/ j4 X" Z1 F
        end- U4 V& ~0 K% ~
    end
    & r$ O8 j$ t% T9 ?alg=pvget(m,'Algorithm');! e7 Q% ], |! c0 G/ O
    maxsize=alg.MaxSize;% K+ k& J" Y$ H
    [A,B,C,D,K,X0]=ssdata(m);
    + _& _# g. s# |8 t1 `# x0 a- ]$ jnx=length(A);
    & c2 U) ?( `- {6 e9 @% g0 h8 Fif ischar(maxsize)
    * i. R9 _1 `1 |( x: u  e7 E5 c    maxsize = idmsize(length(ze{1}),nx);  ^$ _/ Z, C3 ?, I' h1 a
    end
    + a$ K9 p! e1 |6 T0 d7 ~: l0 P5 qAKC=A-K*C;
    ' z3 @8 t1 a: H[ny,nx]=size(C);
    1 a" h6 ?* B4 J% J) Nnu=size(B,2);
    ) L7 _; Q' u! a6 r* z1 R$ `* ~el=zeros(0,ny);
    - r# f$ ~8 }- }( @) i  Pxic = zeros(nx,0);
    & l) }' l; C# p" H8 Q* zNe = length(ze);! m; J9 H, F; e9 R6 m% [0 d2 w
    for kexp = 1:Ne
    / D% G$ m1 |4 M$ W- d" N    z = ze{kexp};7 e2 Q$ v/ O8 Q- A  ^3 k, H9 x
        [Ncap,nz]=size(z);& T5 w0 Z/ Q( t; I& x
        if nu+ny~=nz& Y3 \. ]1 ~  f* L7 k
            error('Ident:general:modelDataDimMismatch',.... u% T( y  a+ Q; n7 e
                'The number of inputs and outputs of the model must match that of the data.')* p. G$ z; O8 `* k
        end* d! J' S* Z' B( T* N9 b
        nz=ny+nu; [Ncap,dum] = size(z); n = nx;
    3 @; Z; \  ?9 M! ^1 L    rowmax = nx+nz; X0 = zeros(nx,1);
    % b% _3 D/ T, Z9 M* _    M=floor(maxsize/rowmax);
    * W7 {7 b% v1 N" E3 P7 L7 b# `    if ny>1 || M<Ncap5 j' L; Y6 P  E
            R=zeros(n,n);Fcap=zeros(n,1);R1=[];
    ) C, [4 v  L) w0 C" x        for kc=1:M:Ncap. C+ w, L" U- ]
                jj=(kc:min(Ncap,kc-1+M));! w+ E$ H/ T+ F; _: y1 S# i
                if jj(length(jj))<Ncap
    8 i4 h$ E9 G4 ~9 \& ^                jjz = [jj,jj(length(jj))+1];; x4 Q2 a7 I9 b- R
                else
    : a/ i: Y* u7 f+ H: m: M                jjz=jj;1 j* L# _2 m/ q! p' H1 J6 r
                end9 y  X) h( u3 T; q
                psitemp=zeros(length(jj),ny);
    ) Z& o) Z4 Y1 m            psi=zeros(ny*length(jj),n);
    5 E& I+ }' M8 Q( z2 S# @' e9 O* U            x=ltitr(AKC,[K B-K*D],z(jjz,:),X0);
    & u  K! z& `  M, G$ U) z            yh=(C*x(1:length(jj),:).').';; f  X! q7 K) ?/ ~) V
                nm=pvget(m,'NoiseVariance');( p6 Q6 @* x4 O! D! S" o* ]
                if isempty(nm) || norm(nm)==0 % To handle models without noise$ i! X! d( {$ h( O: E
                    nm = eye(ny);
    8 }8 K) z/ K/ \* D            end4 N0 Y, o9 N0 @5 I1 I/ t4 j$ a) s5 m
                sqrlam=pinv(sqrtm(nm));
    + ~6 D* e$ a) g. j, }            if ~isempty(D),yh=yh+(D*z(jj,ny+1:ny+nu).').';end
    4 Y# T) ~& E+ {            e=(z(jj,1:ny)-yh)*sqrlam;" x0 R) f8 M' r) W7 T( T
                [nxr,nxc]=size(x);X0=x(nxr,:).';
    ; a) M8 i6 o  K            evec=e(:);6 f+ B$ E1 N7 l3 z: i0 r- f7 \# B
                kl=1;
    - T1 i2 t3 g  }- {2 I            for kx=1:nx
    & t+ T- B& t9 A" S0 X. K1 L                if kc==1
    5 \: c* J, c9 ~  d- |5 V                    x0dum=zeros(nx,1);x0dum(kx,1)=1;
    5 T- f& d' P! p                else
    % z. a# n+ g* S! j                    x0dum=X00(:,kl);9 A; ^3 }7 M+ {& e
                    end
    7 V  V6 Y1 W# A                psix=ltitr(AKC,zeros(nx,1),zeros(length(jjz),1),x0dum);
    3 k6 V5 n( Y8 a# M; ~  v% w% B                [rp,cp]=size(psix);) a' u- u0 {/ q! R+ `- h! F
                    X00(:,kl)=psix(rp,:).';
    5 E7 u1 N6 T; O$ p9 E                psitemp=(C*psix(1:length(jj),:).').'*sqrlam;' b! K" C+ b( c+ |) K
                    psi(:,kl)=psitemp(:);kl=kl+1;) K9 L2 v$ e4 q! `+ a
                end9 V! A5 n& s( J
                if ~isempty(R1)
    5 D. s0 ?9 |& D% b6 y- A$ W                if size(R1,1)<n+1
    4 z+ _+ l" G: S6 Z$ ~                    error('Ident:estimation:predictSmallMaxSize',..." X- T2 i. e0 p0 C
                            'The value of the algorithm property "MaxSize" is too small to estimate the initial states. Either increase MaxSize value in the model or set the "InitialState" option of the "predict" command to a value other than ''Estimate''.')2 L0 ?) o. X! l
                    end( _7 d# V5 e% ^
                    R1=R1(1:n+1,:);
    1 `7 b# _( a0 E& r( h' ^            end
    0 y) e$ f1 K- D0 @5 g- U+ v            H1 = [R1;[psi,evec] ];R1 = triu(qr(H1));0 V9 i# A8 o. s3 S
            end; h6 v4 V9 R$ e" o& o
            try6 s  P; q' I  ]8 v/ S+ r
                xi(:,kexp) = pinv(R1(1:n,1:n))*R1(1:n,n+1);
    2 w9 w+ U  U; D5 g: a( a( R- i  I        catch
    8 Q* ~* F4 z- A            warning('Ident:estimation:X0EstFailed',...
      M  q. z/ s2 r0 r& I. g8 H                'Failed to estimate initial conditions. Check model stability. \nThe "InitialState" option has been set to zero.')
      C) ]5 Q& D3 u* x            xi(:,kexp) = zeros(n,1);
    ' R0 Z1 _+ q" a5 q* s        end
    ! S* l1 x6 c9 y& {) i0 m# u$ w    else! A) Y! P, q- Q8 R. @5 e
            %% First estimate new value of xi' t, e: x: l- N) w) b
            x=ltitr(AKC,[K B-K*D],z);
    8 ^$ R/ p* u& ~3 k' h* \        y0=x*C';
    9 T7 ]* r9 r% \        if ~isempty(D),
    & y/ F2 M8 ?! y! }            y0=y0+(D*z(:,ny+1:ny+nu).').';
    6 c8 Z) b0 o5 S8 D7 {* a2 z        end
    5 b$ G$ u4 {, @6 Y0 T8 |' B; y        psix0=ltitr(AKC.',C.',[1;zeros(Ncap,1)]);3 l2 c* _/ O; u3 P
            psix0=psix0(2:end,:);5 |6 |1 ]* X+ R& E- k) j
            try
    - v$ I  |  I7 v6 X0 E" c            xi(:,kexp) = pinv(psix0)*(z(:,1)-y0);7 p* W) `& \$ ^; u" G7 a4 M! _5 u
            catch& E* |6 H6 J+ P+ D8 d
                warning('Ident:estimation:X0EstFailed',..., \3 z  u! d+ _' M- L
                    'Failed to estimate initial conditions. Check model stability. \nThe "InitialState" option has been set to zero.')
    2 f2 R9 c3 R' f1 G8 g' i3 A/ q            xi(:,kexp) = zeros(n,1);7 I, V" V6 E- A+ P% p1 c5 r' p
            end0 A, l- A# W9 n/ M0 v3 K/ b
        end: k. V7 v: W/ C' W$ V3 u
    end3 E: v7 ^/ O- d
    %%*************************************************************************
    3 U. H0 k  x% ]" o( P) wfunction thpred = polypred(theta,m)
    ) U' u. R; [. [6 k* b: G% Note that cross couplings between output
    5 L2 C  r  a7 \1 ]5 [$ P/ O% channels are ignored in these calculations.
    ) x# p- A& x2 a[ny,nu]= size(theta);9 C3 ^2 e- \4 T  s
    yna = pvget(theta,'OutputName');- \* a2 n% o) x. c. P
    una = pvget(theta,'InputName');7 g6 l2 D+ X7 E( \, L$ P3 z/ |
    yu = pvget(theta,'OutputUnit');! Z' u1 y; Z8 K& r
    uu = pvget(theta,'InputUnit');
    9 Q8 h& C2 F: Q" R3 P! \. cfor ky = 1:ny
    5 ~  ~7 m% ~) J3 ]! z    [a,b,c,d,f] = polydata(theta(ky,:));
    ) J  ^' w' n) ]8 s    if nu>0
    ) K  d0 r0 N& S, y        ff = 1;
    : ?& f+ r" n8 N: w        for ku = 1:nu,
    & Y5 ?, B' j! I            bt = b(ku,:);
    9 L/ S. h2 f7 m, [+ Q& q7 P            for kku = 1:nu, if kku~=ku, bt = conv(bt,f(kku,:)); end, end
    3 U7 K2 n% U( G5 y4 ]2 z. L1 M            bb(ku,:) = bt;  n6 z5 O+ v) K) |9 w6 [( w  @# w
                ff = conv(ff,f(ku,:));
    5 a! g3 W+ x. m& g, v1 c5 I        end, [) ~. r. a4 C# \) t, L
            a = conv(conv(a,ff),d);c=conv(c,ff);; G+ L7 q# C1 G; K# \
        else' G" v  g0 |/ Y& m( u
            a = conv(a,d);
    ; z" }8 C9 \. A  j2 q7 N7 J& y    end3 i7 E# W! C5 ~" j* T( V
        na = length(a); nc = length(c); nn = max(na,nc);
    8 c+ |4 f& p: _7 G$ \6 y    a = [a,zeros(1,nn-na)]; c = [c,zeros(1,nn-nc)];
    ' J7 F! h' n9 s. u( X% t    [f,g] = deconv(conv([1 zeros(1,m-1)],c),a);
    6 \# Q3 F8 Y" W% {& j; M    ng=length(g);3 `) R3 l9 b$ V& J0 C+ O+ r
        if nu>0," S( q$ R2 @9 a* _+ S7 M) y- I
            df=conv(d,f);9 G- ~$ J6 A, b# W8 `  \
            for ku=1:nu
    / @, ^2 f  Q  G0 U! z" w9 g            bf(ku,:)=conv(bb(ku,:),df);
    $ Z5 v0 g2 G9 |! d5 p( H        end
    7 }6 \/ m) Z! e4 W; o5 V4 R1 ^        nbf=length(bf(1,:));nn=max(ng,nbf);4 m/ N, i3 @! u# B# ?: o
            gg=[[g,zeros(1,nn-ng)];[bf,zeros(nu,nn-nbf)]];
    # {+ r8 k( V7 K    else2 Z" l2 Z0 g$ w, y8 r  A
            gg=g;. o1 U3 y& l3 q8 J
        end2 B. \# i7 n0 g0 p
        th1 = idpoly(c,gg);3 a$ K6 g8 z* d6 d  E3 g
        th1 = pvset(th1,'InputName',[yna(ky);una],'OutputName',[yna{ky},'p'],...9 R4 R7 E$ q& Y9 x4 D# S
            'InputUnit',[yu(ky);uu],'OutputUnit',yu(ky),'Ts',pvget(theta,'Ts'),...$ P3 ~" D) a7 i$ B
            'TimeUnit',pvget(theta,'TimeUnit'),'EstimationInfo',pvget(theta,'EstimationInfo'));
    ( v2 _) {% v) y0 Z) t' `2 q. G    thpred{ky} = th1;
    ! M7 U1 u1 l( J; xend
    6 ~* b) l) r$ s4 U%%%************************************************************************
    "真诚赞赏,手留余香"
    还没有人打赏,支持一下
    帖文化:【文明发帖 和谐互动】 社区精神:【创新、交流、互助、共享】
    您需要登录后才可以回帖 登录 | 立即加入

    本版积分规则

    招聘斑竹

    小黑屋|手机版|APP下载(beta)|Archiver|电力研学网 ( 赣ICP备12000811号-1|赣公网安备36040302000210号 )|网站地图

    GMT+8, 2025-4-6 11:50

    Powered by Discuz! X3.5 Licensed

    © 2001-2025 Discuz! Team.

    快速回复 返回顶部 返回列表