Menu
▾
▴
[Octave-cvsupdate] SF.net SVN: octave:[10865] trunk/octave-forge
|
From: <par...@us...> - 2012-08-14 19:59:49
|
Revision: 10865
http://octave.svn.sourceforge.net/octave/?rev=10865&view=rev
Author: paramaniac
Date: 2012-08-14 19:59:40 +0000 (Tue, 14 Aug 2012)
Log Message:
-----------
extra/control-devel ---> main/control
Added Paths:
-----------
trunk/octave-forge/main/control/inst/@iddata/
trunk/octave-forge/main/control/inst/__adjust_iddata__.m
trunk/octave-forge/main/control/inst/__adjust_iddata_tsam__.m
trunk/octave-forge/main/control/inst/__iddata_dim__.m
trunk/octave-forge/main/control/inst/__slicot_identification__.m
trunk/octave-forge/main/control/inst/arx.m
trunk/octave-forge/main/control/inst/moen4.m
trunk/octave-forge/main/control/inst/moesp.m
trunk/octave-forge/main/control/inst/n4sid.m
Removed Paths:
-------------
trunk/octave-forge/extra/control-devel/inst/@iddata/
trunk/octave-forge/extra/control-devel/inst/__adjust_iddata__.m
trunk/octave-forge/extra/control-devel/inst/__adjust_iddata_tsam__.m
trunk/octave-forge/extra/control-devel/inst/__iddata_dim__.m
trunk/octave-forge/extra/control-devel/inst/__slicot_identification__.m
trunk/octave-forge/extra/control-devel/inst/arx.m
trunk/octave-forge/extra/control-devel/inst/moen4.m
trunk/octave-forge/extra/control-devel/inst/moesp.m
trunk/octave-forge/extra/control-devel/inst/n4sid.m
Deleted: trunk/octave-forge/extra/control-devel/inst/__adjust_iddata__.m
===================================================================
--- trunk/octave-forge/extra/control-devel/inst/__adjust_iddata__.m 2012-08-14 19:32:00 UTC (rev 10864)
+++ trunk/octave-forge/extra/control-devel/inst/__adjust_iddata__.m 2012-08-14 19:59:40 UTC (rev 10865)
@@ -1,38 +0,0 @@
-## Copyright (C) 2012 Lukas F. Reichlin
-##
-## This file is part of LTI Syncope.
-##
-## LTI Syncope is free software: you can redistribute it and/or modify
-## it under the terms of the GNU General Public License as published by
-## the Free Software Foundation, either version 3 of the License, or
-## (at your option) any later version.
-##
-## LTI Syncope is distributed in the hope that it will be useful,
-## but WITHOUT ANY WARRANTY; without even the implied warranty of
-## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-## GNU General Public License for more details.
-##
-## You should have received a copy of the GNU General Public License
-## along with LTI Syncope. If not, see <http://www.gnu.org/licenses/>.
-
-## Author: Lukas Reichlin <luk...@gm...>
-## Created: February 2012
-## Version: 0.1
-
-function [y, u] = __adjust_iddata__ (y, u)
-
- if (iscell (y))
- y = reshape (y, [], 1);
- else
- y = {y};
- endif
-
- if (isempty (u))
- u = {}; # avoid [](nx0) and the like
- elseif (iscell (u))
- u = reshape (u, [], 1);
- else
- u = {u};
- endif
-
-endfunction
Deleted: trunk/octave-forge/extra/control-devel/inst/__adjust_iddata_tsam__.m
===================================================================
--- trunk/octave-forge/extra/control-devel/inst/__adjust_iddata_tsam__.m 2012-08-14 19:32:00 UTC (rev 10864)
+++ trunk/octave-forge/extra/control-devel/inst/__adjust_iddata_tsam__.m 2012-08-14 19:59:40 UTC (rev 10865)
@@ -1,52 +0,0 @@
-## Copyright (C) 2012 Lukas F. Reichlin
-##
-## This file is part of LTI Syncope.
-##
-## LTI Syncope is free software: you can redistribute it and/or modify
-## it under the terms of the GNU General Public License as published by
-## the Free Software Foundation, either version 3 of the License, or
-## (at your option) any later version.
-##
-## LTI Syncope is distributed in the hope that it will be useful,
-## but WITHOUT ANY WARRANTY; without even the implied warranty of
-## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-## GNU General Public License for more details.
-##
-## You should have received a copy of the GNU General Public License
-## along with LTI Syncope. If not, see <http://www.gnu.org/licenses/>.
-
-## -*- texinfo -*-
-## Check whether tsam is a e-by-1 cell array of valid sampling times.
-## If not, it tries to convert tsam accordingly.
-## Empty tsam are filled with default value -1.
-
-## Author: Lukas Reichlin <luk...@gm...>
-## Created: February 2012
-## Version: 0.1
-
-function tsam = __adjust_iddata_tsam__ (tsam, e)
-
- if (isempty (tsam))
- tsam = num2cell (-ones (e, 1));
- elseif (iscell (tsam))
- tsam = reshape (tsam, [], 1);
- else
- tsam = {tsam};
- endif
-
- tmp = cellfun (@(x) issample (x, -1), tsam);
-
- if (any (! tmp))
- error ("iddata: invalid sampling time");
- endif
-
- nt = numel (tsam);
-
- if (nt == 1 && e > 1)
- tsam = repmat (tsam, e, 1);
- elseif (nt != e)
- error ("iddata: there are %d experiments, but only %d sampling times", \
- e, nt);
- endif
-
-endfunction
Deleted: trunk/octave-forge/extra/control-devel/inst/__iddata_dim__.m
===================================================================
--- trunk/octave-forge/extra/control-devel/inst/__iddata_dim__.m 2012-08-14 19:32:00 UTC (rev 10864)
+++ trunk/octave-forge/extra/control-devel/inst/__iddata_dim__.m 2012-08-14 19:59:40 UTC (rev 10865)
@@ -1,53 +0,0 @@
-## Author: Lukas Reichlin <luk...@gm...>
-## Created: October 2011
-## Version: 0.1
-
-function [p, m, e] = __iddata_dim__ (y, u)
-
- e = numel (y); # number of experiments
-
- if (isempty (u)) # time series data, outputs only
- [p, m] = cellfun (@__experiment_dim__, y, "uniformoutput", false);
- elseif (e == numel (u)) # outputs and inputs present
- [p, m] = cellfun (@__experiment_dim__, y, u, "uniformoutput", false);
- else
- error ("iddata: require input and output data with matching number of experiments");
- endif
-
- if (e > 1 && ! isequal (p{:}))
- error ("iddata: require identical number of output channels for all experiments");
- endif
-
- if (e > 1 && ! isequal (m{:}))
- error ("iddata: require identical number of input channels for all experiments");
- endif
-
- p = p{1};
- m = m{1};
-
-endfunction
-
-
-function [p, m] = __experiment_dim__ (y, u = [])
-
- if (! is_matrix (y, u))
- error ("iddata: inputs and outputs must be real or complex matrices");
- endif
-
- [ly, p] = size (y);
- [lu, m] = size (u);
-
- if (! isempty (u) && ly != lu)
- error ("iddata: matrices 'y' (%dx%d) and 'u' (%dx%d) must have the same number of samples (rows)", \
- ly, p, lu, m);
- endif
-
- if (ly < p)
- warning ("iddata: more outputs than samples - matrice 'y' should probably be transposed");
- endif
-
- if (lu < m)
- warning ("iddata: more inputs than samples - matrice 'u' should probably be transposed");
- endif
-
-endfunction
Deleted: trunk/octave-forge/extra/control-devel/inst/__slicot_identification__.m
===================================================================
--- trunk/octave-forge/extra/control-devel/inst/__slicot_identification__.m 2012-08-14 19:32:00 UTC (rev 10864)
+++ trunk/octave-forge/extra/control-devel/inst/__slicot_identification__.m 2012-08-14 19:59:40 UTC (rev 10865)
@@ -1,231 +0,0 @@
-## Copyright (C) 2012 Lukas F. Reichlin
-##
-## This file is part of LTI Syncope.
-##
-## LTI Syncope is free software: you can redistribute it and/or modify
-## it under the terms of the GNU General Public License as published by
-## the Free Software Foundation, either version 3 of the License, or
-## (at your option) any later version.
-##
-## LTI Syncope is distributed in the hope that it will be useful,
-## but WITHOUT ANY WARRANTY; without even the implied warranty of
-## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-## GNU General Public License for more details.
-##
-## You should have received a copy of the GNU General Public License
-## along with LTI Syncope. If not, see <http://www.gnu.org/licenses/>.
-
-## -*- texinfo -*-
-## @deftypefn{Function File} {[@var{sys}, @var{x0}], @var{info} =} __slicot_identification__ (@var{method}, @var{dat}, @dots{})
-## Backend for moesp, moen4 and n4sid.
-## @end deftypefn
-
-## Author: Lukas Reichlin <luk...@gm...>
-## Created: May 2012
-## Version: 0.1
-
-function [sys, x0, info] = __slicot_identification__ (method, nout, dat, varargin)
-
- ## determine identification method
- switch (method)
- case "moesp"
- meth = 0;
- case "n4sid"
- meth = 1;
- case "moen4"
- meth = 2;
- otherwise
- error ("ident: invalid method"); # should never happen
- endswitch
-
- if (! isa (dat, "iddata") || ! dat.timedomain)
- error ("%s: first argument must be a time-domain 'iddata' dataset", method);
- endif
-
- if (nargin > 3) # ident (dat, ...)
- if (is_real_scalar (varargin{1})) # ident (dat, n, ...)
- varargin = horzcat (varargin(2:end), {"order"}, varargin(1));
- endif
- if (isstruct (varargin{1})) # ident (dat, opt, ...), ident (dat, n, opt, ...)
- varargin = horzcat (__opt2cell__ (varargin{1}), varargin(2:end));
- endif
- endif
-
- nkv = numel (varargin); # number of keys and values
-
- if (rem (nkv, 2))
- error ("%s: keys and values must come in pairs", method);
- endif
-
- [ns, p, m, e] = size (dat); # dataset dimensions
- tsam = dat.tsam;
-
- ## multi-experiment data requires equal sampling times
- if (e > 1 && ! isequal (tsam{:}))
- error ("%s: require equally sampled experiments", method);
- else
- tsam = tsam{1};
- endif
-
- ## default arguments
- alg = 0;
- conct = 1; # no connection between experiments
- ctrl = 1; # don't confirm order n
- rcond = 0.0;
- tol = 0.0; # -1.0;
- s = [];
- n = [];
- conf = [];
- noise = "n";
-
- ## handle keys and values
- for k = 1 : 2 : nkv
- key = lower (varargin{k});
- val = varargin{k+1};
- switch (key)
- case {"n", "order"}
- if (! issample (val, 0) || val != round (val))
- error ("%s: 'n' must be a positive integer", method);
- endif
- n = val;
- case "s"
- if (! issample (val, 0) || val != round (val))
- error ("%s: 's' must be a positive integer", method);
- endif
- s = val;
- case {"alg", "algorithm"}
- if (strncmpi (val, "c", 1))
- alg = 0; # Cholesky algorithm applied to correlation matrix
- elseif (strncmpi (val, "f", 1))
- alg = 1; # fast QR algorithm
- elseif (strncmpi (val, "q", 1))
- alg = 2; # QR algorithm applied to block Hankel matrices
- else
- error ("%s: invalid algorithm", method);
- endif
- case "tol"
- if (! is_real_scalar (val))
- error ("%s: tolerance 'tol' must be a real scalar", method);
- endif
- tol = val;
- case "rcond"
- if (! is_real_scalar (val))
- error ("%s: 'rcond' must be a real scalar", method);
- endif
- rcond = val;
- case "confirm"
- conf = logical (val);
- case {"input", "inputs"}
- noise = val;
- otherwise
- warning ("%s: invalid property name '%s' ignored", method, key);
- endswitch
- endfor
-
-
- ## handle s/nobr and n
- nsmp = sum (ns); # total number of samples
- nobr = fix ((nsmp+1)/(2*(m+p+1)));
- if (e > 1)
- nobr = min (nobr, fix (min (ns) / 2));
- endif
-
- if (isempty (s) && isempty (n))
- ctrl = 0; # confirm system order estimate
- n = 0;
- elseif (isempty (s))
- s = min (2*n, n+10); # upper bound for n
- nobr = min (nobr, s);
- elseif (isempty (n))
- nobr = __check_s__ (s, nobr, method);
- ctrl = 0; # confirm system order estimate
- n = 0;
- else # s & n non-empty
- nobr = __check_s__ (s, nobr, method);
- if (n >= nobr)
- error ("%s: n=%d, but require n < %d (s)", method, n, nobr);
- endif
- endif
-
- if (! isempty (conf))
- ctrl = ! conf;
- endif
-
- if (nout == 0)
- ## compute singular values
- [sv, nrec] = slib01ad (dat.y, dat.u, nobr, n, meth, alg, conct, ctrl, rcond, tol);
-
- ## there is no 'logbar' function
- svl = log10 (sv);
- base = floor (min (svl));
-
- clf
- bar (svl, "basevalue", base)
- xlim ([0, length(sv)+1])
- yl = ylim;
- ylim ([base, yl(2)])
- title ("Singular Values")
- ylabel ("Logarithm of Singular Values")
- xlabel (sprintf ("Estimated System Order with current Tolerance: %d", nrec))
- grid on
- else
- ## perform system identification
- [a, b, c, d, q, ry, s, k, x0] = slident (dat.y, dat.u, nobr, n, meth, alg, conct, ctrl, rcond, tol);
-
- ## compute noise variance matrix factor L
- ## L L' = Ry, e = L v
- ## v becomes white noise with identity covariance matrix
- l = chol (ry, "lower");
-
- ## assemble model
- [inname, outname] = get (dat, "inname", "outname");
- if (strncmpi (noise, "e", 1)) # add error inputs e, not normalized
- sys = ss (a, [b, k], c, [d, eye(p)], tsam);
- in_u = __labels__ (inname, "u");
- in_e = __labels__ (outname, "y");
- in_e = cellfun (@(x) ["e@", x], in_e, "uniformoutput", false);
- inname = [in_u; in_e];
- elseif (strncmpi (noise, "v", 1)) # add error inputs v, normalized
- sys = ss (a, [b, k*l], c, [d, l], tsam);
- in_u = __labels__ (inname, "u");
- in_v = __labels__ (outname, "y");
- in_v = cellfun (@(x) ["v@", x], in_v, "uniformoutput", false);
- inname = [in_u; in_v];
- elseif (strncmpi (noise, "k", 1)) # Kalman predictor
- sys = ss ([a-k*c], [b-k*d, k], c, [d, zeros(p)], tsam);
- in_u = __labels__ (inname, "u");
- in_y = __labels__ (outname, "y");
- inname = [in_u; in_y];
- else # no error inputs, default
- sys = ss (a, b, c, d, tsam);
- endif
-
- sys = set (sys, "inname", inname, "outname", outname);
-
- ## return x0 as vector for single-experiment data
- ## instead of a cell containing one vector
- if (numel (x0) == 1)
- x0 = x0{1};
- endif
-
- ## assemble info struct
- ## Kalman gain matrix K
- ## state covariance matrix Q
- ## output covariance matrix Ry
- ## state-output cross-covariance matrix S
- ## noise variance matrix factor L
- info = struct ("K", k, "Q", q, "Ry", ry, "S", s, "L", l);
- endif
-
-endfunction
-
-
-function nobr = __check_s__ (s, nobr, method)
-
- if (s <= nobr)
- nobr = s;
- else
- error ("%s: require upper bound s <= %d, but the requested s is %d", method, nobr, s);
- endif
-
-endfunction
Deleted: trunk/octave-forge/extra/control-devel/inst/arx.m
===================================================================
--- trunk/octave-forge/extra/control-devel/inst/arx.m 2012-08-14 19:32:00 UTC (rev 10864)
+++ trunk/octave-forge/extra/control-devel/inst/arx.m 2012-08-14 19:59:40 UTC (rev 10865)
@@ -1,278 +0,0 @@
-## Copyright (C) 2012 Lukas F. Reichlin
-##
-## This file is part of LTI Syncope.
-##
-## LTI Syncope is free software: you can redistribute it and/or modify
-## it under the terms of the GNU General Public License as published by
-## the Free Software Foundation, either version 3 of the License, or
-## (at your option) any later version.
-##
-## LTI Syncope is distributed in the hope that it will be useful,
-## but WITHOUT ANY WARRANTY; without even the implied warranty of
-## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-## GNU General Public License for more details.
-##
-## You should have received a copy of the GNU General Public License
-## along with LTI Syncope. If not, see <http://www.gnu.org/licenses/>.
-
-## -*- texinfo -*-
-## @deftypefn {Function File} {[@var{sys}, @var{x0}] =} arx (@var{dat}, @var{n}, @dots{})
-## @deftypefnx {Function File} {[@var{sys}, @var{x0}] =} arx (@var{dat}, @var{n}, @var{opt}, @dots{})
-## @deftypefnx {Function File} {[@var{sys}, @var{x0}] =} arx (@var{dat}, @var{opt}, @dots{})
-## @deftypefnx {Function File} {[@var{sys}, @var{x0}] =} arx (@var{dat}, @var{'na'}, @var{na}, @var{'nb'}, @var{nb})
-## Estimate ARX model using QR factorization.
-##
-## @strong{Inputs}
-## @table @var
-## @item dat
-## iddata set containing the measurements, i.e. time-domain signals.
-## @item n
-## The desired order of the resulting model @var{sys}.
-## @item @dots{}
-## Optional pairs of keys and values. @code{'key1', value1, 'key2', value2}.
-## @item opt
-## Optional struct with keys as field names.
-## Struct @var{opt} can be created directly or
-## by command @command{options}. @code{opt.key1 = value1, opt.key2 = value2}.
-## @end table
-##
-##
-## @strong{Outputs}
-## @table @var
-## @item sys
-## Discrete-time transfer function model.
-## If the second output argument @var{x0} is returned,
-## @var{sys} becomes a state-space model.
-## @item x0
-## Initial state vector. If @var{dat} is a multi-experiment dataset,
-## @var{x0} becomes a cell vector containing an initial state vector
-## for each experiment.
-## @end table
-##
-##
-## @strong{Option Keys and Values}
-## @table @var
-## @item 'na'
-## Order of the polynomial A(q) and number of poles.
-##
-## @item 'nb'
-## Order of the polynomial B(q)+1 and number of zeros+1.
-## @var{nb} <= @var{na}.
-## @end table
-##
-##
-## @strong{Algorithm}@*
-## Uses the formulae given in [1] on pages 318-319,
-## 'Solving for the LS Estimate by QR Factorization'.
-## For the initial conditions, SLICOT IB01CD is used by courtesy of
-## @uref{http://www.slicot.org, NICONET e.V.}
-##
-## @strong{References}@*
-## [1] Ljung, L. (1999)
-## System Identification - Theory for the User
-## Second Edition
-## Prentice Hall, New Jersey.
-##
-## @end deftypefn
-
-## Author: Lukas Reichlin <luk...@gm...>
-## Created: April 2012
-## Version: 0.1
-
-function [sys, varargout] = arx (dat, varargin)
-
- ## TODO: delays
-
- if (nargin < 2)
- print_usage ();
- endif
-
- if (! isa (dat, "iddata") || ! dat.timedomain)
- error ("arx: first argument must be a time-domain iddata dataset");
- endif
-
- if (is_real_scalar (varargin{1})) # arx (dat, n, ...)
- varargin = horzcat (varargin(2:end), {"na"}, varargin(1), {"nb"}, varargin(1));
- endif
-
- if (isstruct (varargin{1})) # arx (dat, opt, ...), arx (dat, n, opt, ...)
- varargin = horzcat (__opt2cell__ (varargin{1}), varargin(2:end));
- endif
-
- nkv = numel (varargin); # number of keys and values
-
- if (rem (nkv, 2))
- error ("arx: keys and values must come in pairs");
- endif
-
-
- ## p: outputs, m: inputs, ex: experiments
- [~, p, m, ex] = size (dat); # dataset dimensions
-
- ## extract data
- Y = dat.y;
- U = dat.u;
- tsam = dat.tsam;
-
- ## multi-experiment data requires equal sampling times
- if (ex > 1 && ! isequal (tsam{:}))
- error ("arx: require equally sampled experiments");
- else
- tsam = tsam{1};
- endif
-
-
- ## default arguments
- na = [];
- nb = []; % ???
- nk = [];
-
- ## handle keys and values
- for k = 1 : 2 : nkv
- key = lower (varargin{k});
- val = varargin{k+1};
- switch (key)
- ## TODO: proper argument checking
- case "na"
- na = val;
- case "nb"
- nb = val;
- case "nk"
- error ("nk");
- otherwise
- warning ("arx: invalid property name '%s' ignored", key);
- endswitch
- endfor
-
-
- if (is_real_scalar (na, nb))
- na = repmat (na, p, 1); # na(p-by-1)
- nb = repmat (nb, p, m); # nb(p-by-m)
- elseif (! (is_real_vector (na) && is_real_matrix (nb) \
- && rows (na) == p && rows (nb) == p && columns (nb) == m))
- error ("arx: require na(%dx1) instead of (%dx%d) and nb(%dx%d) instead of (%dx%d)", \
- p, rows (na), columns (na), p, m, rows (nb), columns (nb));
- endif
-
- max_nb = max (nb, [], 2); # one maximum for each row/output, max_nb(p-by-1)
- n = max (na, max_nb); # n(p-by-1)
-
- ## create empty cells for numerator and denominator polynomials
- num = cell (p, m+p);
- den = cell (p, m+p);
-
- ## MIMO (p-by-m) models are identified as p MISO (1-by-m) models
- ## For multi-experiment data, minimize the trace of the error
- for i = 1 : p # for every output
- Phi = cell (ex, 1); # one regression matrix per experiment
- for e = 1 : ex # for every experiment
- ## avoid warning: toeplitz: column wins anti-diagonal conflict
- ## therefore set first row element equal to y(1)
- PhiY = toeplitz (Y{e}(1:end-1, i), [Y{e}(1, i); zeros(na(i)-1, 1)]);
- ## create MISO Phi for every experiment
- PhiU = arrayfun (@(x) toeplitz (U{e}(1:end-1, x), [U{e}(1, x); zeros(nb(i,x)-1, 1)]), 1:m, "uniformoutput", false);
- Phi{e} = (horzcat (-PhiY, PhiU{:}))(n(i):end, :);
- endfor
-
- ## compute parameter vector Theta
- Theta = __theta__ (Phi, Y, i, n);
-
- ## extract polynomial matrices A and B from Theta
- ## A is a scalar polynomial for output i, i=1:p
- ## B is polynomial row vector (1-by-m) for output i
- A = [1; Theta(1:na(i))]; # a0 = 1, a1 = Theta(1), an = Theta(n)
- ThetaB = Theta(na(i)+1:end); # all polynomials from B are in one column vector
- B = mat2cell (ThetaB, nb(i,:)); # now separate the polynomials, one for each input
- B = reshape (B, 1, []); # make B a row cell (1-by-m)
- B = cellfun (@(x) [0; x], B, "uniformoutput", false); # b0 = 0 (leading zero required by filt)
-
- ## add error inputs
- Be = repmat ({0}, 1, p); # there are as many error inputs as system outputs (p)
- Be(i) = 1; # inputs m+1:m+p are zero, except m+i which is one
- num(i, :) = [B, Be]; # numerator polynomials for output i, individual for each input
- den(i, :) = repmat ({A}, 1, m+p); # in a row (output i), all inputs have the same denominator polynomial
- endfor
-
- ## A(q) y(t) = B(q) u(t) + e(t)
- ## there is only one A per row
- ## B(z) and A(z) are a Matrix Fraction Description (MFD)
- ## y = A^-1(q) B(q) u(t) + A^-1(q) e(t)
- ## since A(q) is a diagonal polynomial matrix, its inverse is trivial:
- ## the corresponding transfer function has common row denominators.
-
- sys = filt (num, den, tsam); # filt creates a transfer function in z^-1
-
- ## compute initial state vector x0 if requested
- ## this makes only sense for state-space models, therefore convert TF to SS
- if (nargout > 1)
- sys = prescale (ss (sys(:,1:m)));
- x0 = slib01cd (Y, U, sys.a, sys.b, sys.c, sys.d, 0.0);
- ## return x0 as vector for single-experiment data
- ## instead of a cell containing one vector
- if (numel (x0) == 1)
- x0 = x0{1};
- endif
- varargout{1} = x0;
- endif
-
-endfunction
-
-
-function theta = __theta__ (phi, y, i, n)
-
- if (numel (phi) == 1) # single-experiment dataset
- ## use "square-root algorithm"
- A = horzcat (phi{1}, y{1}(n(i)+1:end, i)); # [Phi, Y]
- R0 = triu (qr (A, 0)); # 0 for economy-size R (without zero rows)
- R1 = R0(1:end-1, 1:end-1); # R1 is triangular - can we exploit this in R1\R2?
- R2 = R0(1:end-1, end);
- theta = __ls_svd__ (R1, R2); # R1 \ R2
-
- ## Theta = Phi \ Y(n+1:end, :); # naive formula
- ## theta = __ls_svd__ (phi{1}, y{1}(n(i)+1:end, i));
- else # multi-experiment dataset
- ## TODO: find more sophisticated formula than
- ## Theta = (Phi1' Phi1 + Phi2' Phi2 + ...) \ (Phi1' Y1 + Phi2' Y2 + ...)
-
- ## covariance matrix C = (Phi1' Phi + Phi2' Phi2 + ...)
- tmp = cellfun (@(Phi) Phi.' * Phi, phi, "uniformoutput", false);
- rc = cellfun (@rcond, tmp); # C auch noch testen? QR oder SVD?
- C = plus (tmp{:});
-
- ## PhiTY = (Phi1' Y1 + Phi2' Y2 + ...)
- tmp = cellfun (@(Phi, Y) Phi.' * Y(n(i)+1:end, i), phi, y, "uniformoutput", false);
- PhiTY = plus (tmp{:});
-
- ## pseudoinverse Theta = C \ Phi'Y
- theta = __ls_svd__ (C, PhiTY);
- endif
-
-endfunction
-
-
-function x = __ls_svd__ (A, b)
-
- ## solve the problem Ax=b
- ## x = A\b would also work,
- ## but this way we have better control and warnings
-
- ## solve linear least squares problem by pseudoinverse
- ## the pseudoinverse is computed by singular value decomposition
- ## M = U S V* ---> M+ = V S+ U*
- ## Th = Ph \ Y = Ph+ Y
- ## Th = V S+ U* Y, S+ = 1 ./ diag (S)
-
- [U, S, V] = svd (A, 0); # 0 for "economy size" decomposition
- S = diag (S); # extract main diagonal
- r = sum (S > eps*S(1));
- if (r < length (S))
- warning ("arx: rank-deficient coefficient matrix");
- warning ("sampling time too small");
- warning ("persistence of excitation");
- endif
- V = V(:, 1:r);
- S = S(1:r);
- U = U(:, 1:r);
- x = V * (S .\ (U' * b)); # U' is the conjugate transpose
-
-endfunction
Deleted: trunk/octave-forge/extra/control-devel/inst/moen4.m
===================================================================
--- trunk/octave-forge/extra/control-devel/inst/moen4.m 2012-08-14 19:32:00 UTC (rev 10864)
+++ trunk/octave-forge/extra/control-devel/inst/moen4.m 2012-08-14 19:59:40 UTC (rev 10865)
@@ -1,3034 +0,0 @@
-## Copyright (C) 2012 Lukas F. Reichlin
-##
-## This file is part of LTI Syncope.
-##
-## LTI Syncope is free software: you can redistribute it and/or modify
-## it under the terms of the GNU General Public License as published by
-## the Free Software Foundation, either version 3 of the License, or
-## (at your option) any later version.
-##
-## LTI Syncope is distributed in the hope that it will be useful,
-## but WITHOUT ANY WARRANTY; without even the implied warranty of
-## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-## GNU General Public License for more details.
-##
-## You should have received a copy of the GNU General Public License
-## along with LTI Syncope. If not, see <http://www.gnu.org/licenses/>.
-
-## -*- texinfo -*-
-## @deftypefn {Function File} {[@var{sys}, @var{x0}, @var{info}] =} moen4 (@var{dat}, @dots{})
-## @deftypefnx {Function File} {[@var{sys}, @var{x0}, @var{info}] =} moen4 (@var{dat}, @var{n}, @dots{})
-## @deftypefnx {Function File} {[@var{sys}, @var{x0}, @var{info}] =} moen4 (@var{dat}, @var{opt}, @dots{})
-## @deftypefnx {Function File} {[@var{sys}, @var{x0}, @var{info}] =} moen4 (@var{dat}, @var{n}, @var{opt}, @dots{})
-## Estimate state-space model using combined subspace method:
-## @acronym{MOESP} algorithm for finding the matrices A and C,
-## and @acronym{N4SID} algorithm for finding the matrices B and D.
-## If no output arguments are given, the singular values are
-## plotted on the screen in order to estimate the system order.
-##
-## @strong{Inputs}
-## @table @var
-## @item dat
-## iddata set containing the measurements, i.e. time-domain signals.
-## @item n
-## The desired order of the resulting state-space system @var{sys}.
-## If not specified, @var{n} is chosen automatically according
-## to the singular values and tolerances.
-## @item @dots{}
-## Optional pairs of keys and values. @code{'key1', value1, 'key2', value2}.
-## @item opt
-## Optional struct with keys as field names.
-## Struct @var{opt} can be created directly or
-## by command @command{options}. @code{opt.key1 = value1, opt.key2 = value2}.
-## @end table
-##
-##
-## @strong{Outputs}
-## @table @var
-## @item sys
-## Discrete-time state-space model.
-## @item x0
-## Initial state vector. If @var{dat} is a multi-experiment dataset,
-## @var{x0} becomes a cell vector containing an initial state vector
-## for each experiment.
-## @item info
-## Struct containing additional information.
-## @table @var
-## @item info.K
-## Kalman gain matrix.
-## @item info.Q
-## State covariance matrix.
-## @item info.Ry
-## Output covariance matrix.
-## @item info.S
-## State-output cross-covariance matrix.
-## @item info.L
-## Noise variance matrix factor. LL'=Ry.
-## @end table
-## @end table
-##
-##
-##
-## @strong{Option Keys and Values}
-## @table @var
-## @item 'n'
-## The desired order of the resulting state-space system @var{sys}.
-## @var{s} > @var{n} > 0.
-##
-## @item 's'
-## The number of block rows @var{s} in the input and output
-## block Hankel matrices to be processed. @var{s} > 0.
-## In the MOESP theory, @var{s} should be larger than @var{n},
-## the estimated dimension of state vector.
-##
-## @item 'alg', 'algorithm'
-## Specifies the algorithm for computing the triangular
-## factor R, as follows:
-## @table @var
-## @item 'C'
-## Cholesky algorithm applied to the correlation
-## matrix of the input-output data. Default method.
-## @item 'F'
-## Fast QR algorithm.
-## @item 'Q'
-## QR algorithm applied to the concatenated block
-## Hankel matrices.
-## @end table
-##
-## @item 'tol'
-## Absolute tolerance used for determining an estimate of
-## the system order. If @var{tol} >= 0, the estimate is
-## indicated by the index of the last singular value greater
-## than or equal to @var{tol}. (Singular values less than @var{tol}
-## are considered as zero.) When @var{tol} = 0, an internally
-## computed default value, @var{tol} = @var{s}*@var{eps}*SV(1), is used,
-## where SV(1) is the maximal singular value, and @var{eps} is
-## the relative machine precision.
-## When @var{tol} < 0, the estimate is indicated by the
-## index of the singular value that has the largest
-## logarithmic gap to its successor. Default value is 0.
-##
-## @item 'rcond'
-## The tolerance to be used for estimating the rank of
-## matrices. If the user sets @var{rcond} > 0, the given value
-## of @var{rcond} is used as a lower bound for the reciprocal
-## condition number; an m-by-n matrix whose estimated
-## condition number is less than 1/@var{rcond} is considered to
-## be of full rank. If the user sets @var{rcond} <= 0, then an
-## implicitly computed, default tolerance, defined by
-## @var{rcond} = m*n*@var{eps}, is used instead, where @var{eps} is the
-## relative machine precision. Default value is 0.
-##
-## @item 'confirm'
-## Specifies whether or not the user's confirmation of the
-## system order estimate is desired, as follows:
-## @table @var
-## @item true
-## User's confirmation.
-## @item false
-## No confirmation. Default value.
-## @end table
-##
-## @item 'input'
-## The desired type of noise input channels.
-## @table @var
-## @item 'n'
-## No error inputs. Default value.
-## @iftex
-## @tex
-## $$ x_{k+1} = A x_k + B u_k $$
-## $$ y_k = C x_k + D u_k $$
-## @end tex
-## @end iftex
-## @ifnottex
-## @example
-## x[k+1] = A x[k] + B u[k]
-## y[k] = C x[k] + D u[k]
-## @end example
-## @end ifnottex
-##
-## @item 'e'
-## Return @var{sys} as a (p-by-m+p) state-space model with
-## both measured input channels u and noise channels e
-## with covariance matrix @var{Ry}.
-## @iftex
-## @tex
-## $$ x_{k+1} = A x_k + B u_k + K e_k $$
-## $$ y_k = C x_k + D u_k + e_k $$
-## @end tex
-## @end iftex
-## @ifnottex
-## @example
-## x[k+1] = A x[k] + B u[k] + K e[k]
-## y[k] = C x[k] + D u[k] + e[k]
-## @end example
-## @end ifnottex
-##
-## @item 'v'
-## Return @var{sys} as a (p-by-m+p) state-space model with
-## both measured input channels u and white noise channels v
-## with identity covariance matrix.
-## @iftex
-## @tex
-## $$ x_{k+1} = A x_k + B u_k + K L v_k $$
-## $$ y_k = C x_k + D u_k + L v_k $$
-## $$ e = L v, \\ L L^T = R_y $$
-## @end tex
-## @end iftex
-## @ifnottex
-## @example
-## x[k+1] = A x[k] + B u[k] + K L v[k]
-## y[k] = C x[k] + D u[k] + L v[k]
-## e = L v, L L' = Ry
-## @end example
-## @end ifnottex
-##
-## @item 'k'
-## Return @var{sys} as a Kalman predictor for simulation.
-## @iftex
-## @tex
-## $$ \\widehat{x}_{k+1} = A \\widehat{x}_k + B u_k + K (y_k - \\widehat{y}_k) $$
-## $$ \\widehat{y}_k = C \\widehat{x}_k + D u_k $$
-## @end tex
-## @end iftex
-## @ifnottex
-## @example
-## ^ ^ ^
-## x[k+1] = A x[k] + B u[k] + K(y[k] - y[k])
-## ^ ^
-## y[k] = C x[k] + D u[k]
-## @end example
-## @end ifnottex
-##
-## @iftex
-## @tex
-## $$ \\widehat{x}_{k+1} = (A-KC) \\widehat{x}_k + (B-KD) u_k + K y_k $$
-## $$ \\widehat{y}_k = C \\widehat{x}_k + D u_k + 0 y_k $$
-## @end tex
-## @end iftex
-## @ifnottex
-## @example
-## ^ ^
-## x[k+1] = (A-KC) x[k] + (B-KD) u[k] + K y[k]
-## ^ ^
-## y[k] = C x[k] + D u[k] + 0 y[k]
-## @end example
-## @end ifnottex
-## @end table
-## @end table
-##
-##
-## @strong{Algorithm}@*
-## Uses SLICOT IB01AD, IB01BD and IB01CD by courtesy of
-## @uref{http://www.slicot.org, NICONET e.V.}
-##
-## @end deftypefn
-
-## Author: Lukas Reichlin <luk...@gm...>
-## Created: May 2012
-## Version: 0.1
-
-function [sys, x0, info] = moen4 (varargin)
-
- if (nargin == 0)
- print_usage ();
- endif
-
- if (nargout == 0)
- __slicot_identification__ ("moen4", nargout, varargin{:});
- else
- [sys, x0, info] = __slicot_identification__ ("moen4", nargout, varargin{:});
- endif
-
-endfunction
-
-
-%!shared SYS, X0, INFO, Ae, Be, Ce, De, Ke, Qe, Rye, Se, X0e
-%!
-%! Y = [ 4.7661 5.5451 5.8503 5.3766 4.8833 5.4865 3.5378 5.3155 6.0530 4.3729
-%! 4.7637 5.1886 5.9236 5.6818 4.8858 5.1495 3.5549 5.5329 6.0799 4.7417
-%! 4.8394 4.8833 5.9212 5.8235 4.8931 4.8442 3.4938 5.4450 6.1287 5.0884
-%! 5.0030 4.6000 5.9773 5.9529 4.7148 4.5414 3.4474 5.3961 6.0799 5.1861
-%! 5.0176 4.2704 5.7405 6.0628 4.4511 4.2679 3.4401 5.2740 6.1678 5.0372
-%! 5.0567 4.0384 5.3888 6.0897 4.2337 4.0604 3.4083 5.0274 6.1947 4.7856
-%! 5.1544 3.8381 5.0005 6.0750 4.0433 3.9602 3.4108 4.7441 6.2362 4.5634
-%! 5.3619 3.7112 4.8491 6.0262 3.8650 3.7893 3.4523 4.6684 6.0530 4.5341
-%! 5.4254 3.5915 4.9444 5.9944 3.7576 3.6428 3.6818 4.6513 5.6525 4.7050
-%! 5.5695 3.5353 5.1739 6.0775 3.6696 3.5256 4.0604 4.5146 5.2740 4.7417
-%! 5.6818 3.4865 5.3693 5.8577 3.5939 3.4987 4.4413 4.2679 4.8589 4.6489
-%! 5.7429 3.4767 5.4474 5.7014 3.5475 3.4547 4.8540 4.2606 4.5341 4.4315
-%! 5.8039 3.4254 5.6037 5.7307 3.5060 3.4083 5.1544 4.2630 4.4560 4.2386
-%! 5.9187 3.3815 5.7307 5.7844 3.4547 3.3790 5.4254 4.1898 4.6196 4.0652
-%! 5.8210 3.3693 5.8503 5.8235 3.3986 3.3766 5.5964 4.2777 4.8662 3.9431
-%! 5.4474 3.3644 5.9798 5.8943 3.3619 3.3619 5.5866 4.6000 5.1177 3.8113
-%! 5.0616 3.3473 5.9920 5.7624 3.3400 3.3595 5.3546 4.9322 5.1666 3.6916
-%! 4.6293 3.3815 6.0848 5.4157 3.3742 3.3693 5.0274 5.2838 5.0567 3.6525
-%! 4.2679 3.4206 5.9407 4.9615 3.5207 3.3986 4.8638 5.5280 5.0030 3.8259
-%! 4.0115 3.4132 5.8039 4.5952 3.7136 3.5793 4.7612 5.7405 5.0982 4.2240
-%! 3.8503 3.4523 5.7917 4.3314 3.7576 3.9480 4.5707 5.8748 5.3253 4.4242
-%! 3.7112 3.6355 5.6037 4.2972 3.7795 4.4120 4.3681 5.9554 5.5671 4.4291
-%! 3.5695 4.0384 5.2643 4.5829 3.6965 4.5854 4.3974 5.9920 5.4670 4.3192
-%! 3.5182 4.3754 4.9468 4.8613 3.7771 4.5146 4.5732 5.8455 5.2521 4.1385
-%! 3.6525 4.7270 4.6196 5.1739 3.8870 4.3436 4.8418 5.5280 4.9468 3.9651
-%! 3.8186 5.0567 4.5146 5.1666 3.9041 4.1556 5.2032 5.0616 4.8809 3.8870
-%! 3.8626 5.2985 4.4340 4.9199 3.8503 3.9847 5.4523 4.7344 4.9810 3.8015
-%! 4.0115 5.5329 4.2850 4.6074 3.9651 4.0433 5.6525 4.5341 5.2252 3.7014
-%! 4.3534 5.4670 4.1214 4.3705 4.2826 4.3070 5.8552 4.5341 5.4596 3.6403
-%! 4.7050 5.1959 3.9456 4.1825 4.5219 4.4218 5.9065 4.6977 5.7234 3.7673
-%! 5.0836 4.8858 3.9847 4.0384 4.7148 4.3534 5.9529 4.7441 5.7917 4.1507
-%! 5.3449 4.7637 4.2191 4.1458 4.9712 4.2240 5.8284 4.6196 5.9065 4.6489
-%! 5.2740 4.8760 4.5463 4.4315 5.2203 4.0530 5.7917 4.6440 5.9920 4.9908
-%! 5.1275 5.0420 4.8735 4.5561 5.5329 3.9407 5.7991 4.8320 5.8357 5.0884
-%! 4.7612 5.2838 5.1544 4.4804 5.6525 3.8381 5.8137 5.1324 5.5280 5.0225
-%! 4.4511 5.4914 5.3888 4.3754 5.7820 3.7307 5.8772 5.4108 5.1422 4.7832
-%! 4.2215 5.5964 5.6135 4.3705 5.9554 3.6525 5.9554 5.6257 4.7759 4.6855
-%! 4.0457 5.6721 5.8357 4.5585 6.0359 3.6110 5.7820 5.6037 4.4902 4.6660
-%! 3.8748 5.7722 5.8845 4.8589 6.1190 3.5646 5.5182 5.3155 4.2362 4.7075
-%! 3.7307 5.8308 5.9554 4.8955 6.1336 3.4963 5.1275 4.9615 4.0237 4.9126
-%! 3.6623 5.9334 5.7624 4.7417 6.1532 3.4621 4.7637 4.6196 3.8870 5.1959
-%! 3.5768 5.8992 5.4596 4.7441 6.1922 3.4547 4.4926 4.3583 3.7527 5.4157
-%! 3.5427 5.9358 5.0616 4.8760 6.1434 3.4254 4.2337 4.1556 3.6818 5.6232
-%! 3.4792 5.8943 4.7075 5.1055 6.1678 3.3790 4.0115 4.0335 3.8064 5.7405
-%! 3.4547 5.9187 4.4584 5.2398 5.9920 3.4328 3.8552 3.8870 4.1458 5.8992
-%! 3.3595 5.9944 4.2679 5.5182 5.6525 3.6232 3.6916 3.7722 4.6000 5.9285
-%! 3.2985 5.9578 4.0530 5.6525 5.4596 3.9749 3.6355 3.6403 5.0030 6.0506
-%! 3.2252 6.0311 3.9431 5.7234 5.4376 4.3803 3.8186 3.5329 5.3033 6.1532
-%! 3.2008 6.0628 3.8259 5.8552 5.3400 4.7148 4.1556 3.4352 5.5524 5.9651
-%! 3.2252 6.0408 3.9676 5.9627 5.0982 5.0738 4.5903 3.4279 5.6159 5.5866
-%! 3.2276 6.0970 4.2801 5.9847 4.7856 5.3693 4.9883 3.4230 5.5231 5.3815
-%! 3.2740 6.1239 4.4804 5.9847 4.4926 5.6037 5.0762 3.3986 5.6110 5.3717
-%! 3.4572 6.1629 4.4926 6.0555 4.2362 5.7453 4.9077 3.6037 5.7136 5.4865
-%! 3.8674 6.0408 4.3900 6.0628 4.0677 5.6525 4.6489 4.0237 5.8455 5.5671
-%! 4.3217 5.8455 4.1971 6.0555 3.9334 5.4010 4.3778 4.4511 5.8992 5.8210
-%! 4.4926 5.7722 4.1116 6.0701 3.8235 5.0152 4.2166 4.7930 5.9944 5.9138
-%! 4.4315 5.7991 3.9822 5.7844 3.7307 4.7099 4.2875 4.9029 6.0921 5.9944
-%! 4.2435 5.9236 3.8674 5.4401 3.6110 4.4169 4.5903 4.7808 6.0921 6.0115
-%! 4.0506 5.9285 3.7673 5.0567 3.5646 4.2362 4.8467 4.5903 6.1434 5.9993
-%! 3.8577 6.0018 3.8723 4.9419 3.5500 4.2362 5.1397 4.3363 6.1532 6.0188
-%! 3.7307 6.0018 4.2362 5.0103 3.5573 4.2484 5.3888 4.1458 6.2337 5.8210
-%! 3.7917 6.0604 4.6635 5.1348 3.5134 4.2215 5.6892 4.2166 6.1873 5.7282
-%! 3.9212 5.8821 4.9712 5.3131 3.5158 4.2972 5.8845 4.4340 6.0140 5.7405
-%! 3.9554 5.5109 5.0665 5.4792 3.6941 4.5903 6.0433 4.7148 5.8357 5.7649
-%! 3.8479 5.3229 4.9029 5.6232 4.0726 4.8931 6.1703 5.0982 5.7746 5.8821
-%! 3.7258 5.3717 4.6757 5.5622 4.4804 5.1348 6.2118 5.3595 5.6867 5.9260
-%! 3.6110 5.4547 4.3925 5.3302 4.7050 5.4279 6.2508 5.5695 5.5378 5.7502
-%! 3.7160 5.4376 4.0994 5.0103 4.6123 5.3790 6.2093 5.7722 5.3278 5.4157
-%! 4.0921 5.1593 4.1141 4.6660 4.3851 5.3644 6.0140 5.9212 5.0543 4.9956
-%! 4.4804 4.9029 4.3265 4.4145 4.2020 5.4523 5.7014 6.0555 4.7002 4.8613
-%! 4.8149 4.5878 4.6440 4.2020 4.0262 5.5671 5.4694 5.9627 4.3949 4.9029
-%! 5.0543 4.5024 4.9712 4.0482 3.9041 5.6721 5.4792 5.6428 4.1800 5.1031
-%! 5.3033 4.5952 5.1593 4.0799 3.7746 5.7698 5.5573 5.4352 4.0433 5.3644
-%! 5.4865 4.8247 5.3888 4.1898 3.6916 5.8308 5.7282 5.3888 3.8772 5.5964
-%! 5.6721 5.0640 5.5768 4.1312 3.8455 5.9236 5.8821 5.5378 3.7527 5.7527
-%! 5.7795 5.2716 5.6525 4.0042 4.2020 5.9651 5.9847 5.6818 3.7282 5.8455
-%! 5.7991 5.4670 5.8039 3.9163 4.5854 6.0579 5.9016 5.7014 3.8699 5.9285
-%! 5.6648 5.6159 5.9138 3.9602 4.9029 6.0506 5.5817 5.6159 4.2069 6.0066
-%! 5.2911 5.5280 5.8870 4.1996 5.2569 6.0726 5.3717 5.6672 4.3558 5.8406
-%! 4.8809 5.2545 5.7991 4.6245 5.5109 6.1116 5.4181 5.7405 4.4267 5.5182
-%! 4.5585 4.8833 5.7307 4.8833 5.6403 6.0701 5.5109 5.8039 4.4535 5.1739
-%! 4.1849 4.5170 5.7624 5.1373 5.8430 5.8967 5.6672 5.8821 4.5219 4.7392
-%! 3.8894 4.1971 5.8137 5.3790 5.9749 5.7551 5.7917 5.9505 4.3925 4.4584
-%! 3.7087 4.0018 5.8210 5.6232 5.9358 5.7185 5.6989 6.0726 4.1556 4.4267
-%! 3.6232 3.8064 5.9285 5.7624 5.8210 5.8210 5.4840 6.1483 3.9651 4.6025
-%! 3.5695 3.9041 6.0140 5.8333 5.5280 6.0018 5.1544 6.1165 3.8772 4.8223
-%! 3.7185 3.9236 5.7649 5.6867 5.1715 6.0018 4.9810 6.1776 3.9700 5.1837
-%! 4.0335 3.8699 5.4132 5.3668 4.8101 5.9016 5.0616 6.2020 4.2582 5.4303
-%! 4.4120 3.8064 5.0982 5.2252 4.4535 5.5573 5.1959 6.2069 4.4218 5.6525
-%! 4.6293 3.7209 4.6782 5.2398 4.3803 5.1739 5.3595 5.9920 4.3363 5.8210
-%! 4.5585 3.8186 4.3729 5.3546 4.5659 4.8003 5.6159 5.5646 4.2997 5.7063
-%! 4.3949 4.1409 4.3925 5.5085 4.8052 4.4315 5.7624 5.1788 4.3925 5.3693
-%! 4.1800 4.5292 4.5903 5.5964 5.1251 4.1947 5.8577 4.9981 4.6757 5.0274
-%! 4.1971 4.8052 4.9199 5.7527 5.3546 4.0066 5.9480 5.0518 4.7612 4.7050
-%! 4.4315 5.0860 5.0176 5.8748 5.5891 3.8503 5.8357 5.2325 4.6587 4.4145
-%! 4.7148 5.3400 4.8589 5.9065 5.7649 3.7478 5.7063 5.4840 4.4902 4.1458
-%! 4.9615 5.5329 4.6757 5.8943 5.9236 3.6428 5.4987 5.6867 4.3070 3.9651
-%! 5.3009 5.5768 4.6196 5.7429 5.9407 3.5915 5.1886 5.8992 4.1263 4.0335
-%! 5.5671 5.6672 4.8345 5.4474 5.8577 3.5695 5.1177 5.8699 3.9724 4.3729
-%! 5.6818 5.7917 5.0909 5.0250 5.6941 3.5280 5.1910 5.9773 4.0775 4.6831 ](:);
-%!
-%!
-%!
-%! U = [ 6.4100 3.4100 6.4100 6.4100 3.4100 3.4100 3.4100 6.4100 6.4100 3.4100
-%! 3.4100 3.4100 3.4100 6.4100 3.4100 3.4100 3.4100 3.4100 6.4100 3.4100
-%! 6.4100 3.4100 3.4100 6.4100 3.4100 3.4100 3.4100 3.4100 6.4100 3.4100
-%! 6.4100 3.4100 3.4100 6.4100 3.4100 3.4100 3.4100 3.4100 6.4100 3.4100
-%! 6.4100 3.4100 6.4100 6.4100 3.4100 3.4100 3.4100 6.4100 3.4100 6.4100
-%! 6.4100 3.4100 6.4100 6.4100 3.4100 3.4100 6.4100 3.4100 3.4100 6.4100
-%! 6.4100 3.4100 6.4100 6.4100 3.4100 3.4100 6.4100 3.4100 3.4100 3.4100
-%! 6.4100 3.4100 6.4100 3.4100 3.4100 3.4100 6.4100 3.4100 3.4100 3.4100
-%! 6.4100 3.4100 6.4100 6.4100 3.4100 3.4100 6.4100 6.4100 3.4100 3.4100
-%! 6.4100 3.4100 6.4100 6.4100 3.4100 3.4100 6.4100 3.4100 6.4100 3.4100
-%! 6.4100 3.4100 6.4100 6.4100 3.4100 3.4100 6.4100 3.4100 6.4100 3.4100
-%! 3.4100 3.4100 6.4100 6.4100 3.4100 3.4100 6.4100 6.4100 6.4100 3.4100
-%! 3.4100 3.4100 6.4100 6.4100 3.4100 3.4100 3.4100 6.4100 6.4100 3.4100
-%! 3.4100 3.4100 6.4100 3.4100 3.4100 3.4100 3.4100 6.4100 3.4100 3.4100
-%! 3.4100 3.4100 6.4100 3.4100 3.4100 3.4100 3.4100 6.4100 3.4100 3.4100
-%! 3.4100 3.4100 3.4100 3.4100 6.4100 3.4100 6.4100 6.4100 6.4100 6.4100
-%! 3.4100 3.4100 6.4100 3.4100 3.4100 6.4100 3.4100 6.4100 6.4100 6.4100
-%! 3.4100 3.4100 6.4100 3.4100 3.4100 6.4100 3.4100 6.4100 6.4100 3.4100
-%! 3.4100 6.4100 3.4100 6.4100 3.4100 6.4100 3.4100 6.4100 6.4100 3.4100
-%! 3.4100 6.4100 3.4100 6.4100 3.4100 3.4100 6.4100 6.4100 3.4100 3.4100
-%! 3.4100 6.4100 3.4100 6.4100 6.4100 3.4100 6.4100 3.4100 3.4100 3.4100
-%! 6.4100 6.4100 3.4100 6.4100 3.4100 3.4100 6.4100 3.4100 3.4100 3.4100
-%! 3.4100 6.4100 6.4100 3.4100 3.4100 3.4100 6.4100 3.4100 6.4100 3.4100
-%! 3.4100 6.4100 3.4100 3.4100 3.4100 3.4100 6.4100 3.4100 6.4100 3.4100
-%! 6.4100 6.4100 3.4100 3.4100 6.4100 6.4100 6.4100 3.4100 6.4100 3.4100
-%! 6.4100 3.4100 3.4100 3.4100 6.4100 6.4100 6.4100 6.4100 6.4100 3.4100
-%! 6.4100 3.4100 3.4100 3.4100 3.4100 3.4100 6.4100 6.4100 6.4100 6.4100
-%! 6.4100 3.4100 6.4100 3.4100 6.4100 3.4100 6.4100 3.4100 6.4100 6.4100
-%! 6.4100 6.4100 6.4100 6.4100 6.4100 3.4100 3.4100 3.4100 6.4100 6.4100
-%! 3.4100 6.4100 6.4100 6.4100 6.4100 3.4100 6.4100 6.4100 6.4100 6.4100
-%! 3.4100 6.4100 6.4100 3.4100 6.4100 3.4100 6.4100 6.4100 3.4100 3.4100
-%! 3.4100 6.4100 6.4100 3.4100 6.4100 3.4100 6.4100 6.4100 3.4100 3.4100
-%! 3.4100 6.4100 6.4100 3.4100 6.4100 3.4100 6.4100 6.4100 3.4100 3.4100
-%! 3.4100 6.4100 6.4100 6.4100 6.4100 3.4100 6.4100 6.4100 3.4100 6.4100
-%! 3.4100 6.4100 6.4100 6.4100 6.4100 3.4100 3.4100 3.4100 3.4100 3.4100
-%! 3.4100 6.4100 6.4100 6.4100 6.4100 3.4100 3.4100 3.4100 3.4100 6.4100
-%! 3.4100 6.4100 6.4100 3.4100 6.4100 3.4100 3.4100 3.4100 3.4100 6.4100
-%! 3.4100 6.4100 3.4100 3.4100 6.4100 3.4100 3.4100 3.4100 3.4100 6.4100
-%! 3.4100 6.4100 3.4100 6.4100 6.4100 3.4100 3.4100 3.4100 3.4100 6.4100
-%! 3.4100 6.4100 3.4100 6.4100 6.4100 3.4100 3.4100 3.4100 3.4100 6.4100
-%! 3.4100 6.4100 3.4100 6.4100 6.4100 3.4100 3.4100 3.4100 6.4100 6.4100
-%! 3.4100 6.4100 3.4100 6.4100 3.4100 3.4100 3.4100 3.4100 6.4100 6.4100
-%! 3.4100 6.4100 3.4100 6.4100 3.4100 6.4100 3.4100 3.4100 6.4100 6.4100
-%! 3.4100 6.4100 3.4100 6.4100 6.4100 6.4100 3.4100 3.4100 6.4100 6.4100
-%! 3.4100 6.4100 3.4100 6.4100 6.4100 6.4100 6.4100 3.4100 6.4100 6.4100
-%! 3.4100 6.4100 3.4100 6.4100 3.4100 6.4100 6.4100 3.4100 6.4100 3.4100
-%! 3.4100 6.4100 6.4100 6.4100 3.4100 6.4100 6.4100 3.4100 3.4100 3.4100
-%! 3.4100 6.4100 6.4100 6.4100 3.4100 6.4100 6.4100 3.4100 6.4100 6.4100
-%! 3.4100 6.4100 3.4100 6.4100 3.4100 6.4100 3.4100 3.4100 6.4100 6.4100
-%! 6.4100 6.4100 3.4100 6.4100 3.4100 6.4100 3.4100 6.4100 6.4100 6.4100
-%! 6.4100 3.4100 3.4100 6.4100 3.4100 3.4100 3.4100 6.4100 6.4100 6.4100
-%! 6.4100 6.4100 3.4100 6.4100 3.4100 3.4100 3.4100 6.4100 6.4100 6.4100
-%! 3.4100 6.4100 3.4100 6.4100 3.4100 3.4100 3.4100 6.4100 6.4100 6.4100
-%! 3.4100 6.4100 3.4100 3.4100 3.4100 3.4100 6.4100 3.4100 6.4100 6.4100
-%! 3.4100 6.4100 3.4100 3.4100 3.4100 3.4100 6.4100 3.4100 6.4100 6.4100
-%! 3.4100 6.4100 3.4100 3.4100 3.4100 3.4100 6.4100 3.4100 6.4100 6.4100
-%! 3.4100 6.4100 6.4100 6.4100 3.4100 6.4100 6.4100 3.4100 6.4100 6.4100
-%! 3.4100 6.4100 6.4100 6.4100 3.4100 3.4100 6.4100 3.4100 6.4100 3.4100
-%! 6.4100 6.4100 6.4100 6.4100 3.4100 3.4100 6.4100 6.4100 6.4100 6.4100
-%! 3.4100 3.4100 6.4100 6.4100 3.4100 6.4100 6.4100 6.4100 3.4100 6.4100
-%! 3.4100 3.4100 3.4100 6.4100 6.4100 6.4100 6.4100 6.4100 6.4100 6.4100
-%! 3.4100 6.4100 3.4100 6.4100 6.4100 6.4100 6.4100 6.4100 6.4100 6.4100
-%! 3.4100 6.4100 3.4100 3.4100 6.4100 6.4100 6.4100 6.4100 3.4100 6.4100
-%! 3.4100 6.4100 3.4100 3.4100 3.4100 6.4100 6.4100 6.4100 6.4100 3.4100
-%! 6.4100 3.4100 3.4100 3.4100 3.4100 3.4100 6.4100 6.4100 3.4100 3.4100
-%! 6.4100 3.4100 6.4100 3.4100 3.4100 6.4100 3.4100 6.4100 3.4100 3.4100
-%! 6.4100 3.4100 6.4100 3.4100 3.4100 6.4100 3.4100 6.4100 3.4100 6.4100
-%! 6.4100 3.4100 6.4100 3.4100 3.4100 6.4100 6.4100 3.4100 3.4100 6.4100
-%! 6.4100 6.4100 6.4100 3.4100 3.4100 6.4100 6.4100 3.4100 3.4100 6.4100
-%! 6.4100 6.4100 6.4100 6.4100 3.4100 6.4100 6.4100 6.4100 3.4100 6.4100
-%! 6.4100 6.4100 6.4100 3.4100 3.4100 6.4100 6.4100 6.4100 3.4100 6.4100
-%! 6.4100 6.4100 6.4100 3.4100 6.4100 6.4100 6.4100 6.4100 3.4100 6.4100
-%! 6.4100 6.4100 6.4100 3.4100 6.4100 6.4100 6.4100 6.4100 3.4100 6.4100
-%! 6.4100 6.4100 6.4100 3.4100 6.4100 6.4100 3.4100 3.4100 6.4100 6.4100
-%! 3.4100 6.4100 6.4100 6.4100 6.4100 6.4100 3.4100 6.4100 6.4100 6.4100
-%! 3.4100 3.4100 6.4100 6.4100 6.4100 6.4100 6.4100 6.4100 3.4100 3.4100
-%! 3.4100 3.4100 3.4100 6.4100 6.4100 6.4100 6.4100 6.4100 3.4100 3.4100
-%! 3.4100 3.4100 6.4100 6.4100 6.4100 6.4100 6.4100 6.4100 6.4100 3.4100
-%! 3.4100 3.4100 6.4100 6.4100 6.4100 3.4100 6.4100 6.4100 3.4100 3.4100
-%! 3.4100 3.4100 6.4100 6.4100 6.4100 6.4100 6.4100 6.4100 3.4100 3.4100
-%! 3.4100 3.4100 6.4100 6.4100 6.4100 6.4100 3.4100 6.4100 3.4100 6.4100
-%! 3.4100 3.4100 6.4100 6.4100 3.4100 6.4100 3.4100 6.4100 3.4100 6.4100
-%! 3.4100 6.4100 6.4100 6.4100 3.4100 6.4100 3.4100 6.4100 3.4100 6.4100
-%! 6.4100 3.4100 3.4100 3.4100 3.4100 6.4100 6.4100 6.4100 6.4100 6.4100
-%! 6.4100 3.4100 3.4100 3.4100 3.4100 3.4100 6.4100 6.4100 6.4100 6.4100
-%! 6.4100 3.4100 3.4100 6.4100 3.4100 3.4100 6.4100 6.4100 3.4100 6.4100
-%! 3.4100 3.4100 3.4100 6.4100 6.4100 3.4100 6.4100 3.4100 3.4100 6.4100
-%! 3.4100 6.4100 3.4100 6.4100 6.4100 3.4100 6.4100 3.4100 3.4100 3.4100
-%! 3.4100 6.4100 6.4100 6.4100 6.4100 3.4100 6.4100 3.4100 6.4100 3.4100
-%! 3.4100 6.4100 6.4100 6.4100 6.4100 3.4100 6.4100 6.4100 6.4100 3.4100
-%! 6.4100 6.4100 6.4100 6.4100 6.4100 3.4100 6.4100 6.4100 3.4100 3.4100
-%! 6.4100 6.4100 3.4100 6.4100 6.4100 3.4100 3.4100 6.4100 3.4100 3.4100
-%! 6.4100 6.4100 3.4100 6.4100 6.4100 3.4100 6.4100 6.4100 3.4100 3.4100
-%! 6.4100 6.4100 3.4100 6.4100 6.4100 3.4100 3.4100 6.4100 3.4100 3.4100
-%! 6.4100 6.4100 6.4100 3.4100 6.4100 3.4100 3.4100 6.4100 3.4100 6.4100
-%! 6.4100 6.4100 6.4100 3.4100 3.4100 3.4100 6.4100 6.4100 3.4100 6.4100
-%! 6.4100 6.4100 6.4100 3.4100 6.4100 3.4100 6.4100 6.4100 6.4100 6.4100
-%! 3.4100 6.4100 6.4100 6.4100 3.4100 3.4100 6.4100 6.4100 6.4100 6.4100
-%! 3.4100 6.4100 6.4100 6.4100 3.4100 3.4100 6.4100 6.4100 6.4100 6.4100
-%! 3.4100 6.4100 6.4100 3.4100 3.4100 3.4100 3.4100 6.4100 6.4100 6.4100 ](:);
-%!
-%!
-%! DAT = iddata (Y, U);
-%!
-%! [SYS, X0, INFO] = moen4 (DAT, "s", 15, "rcond", 0.0, "tol", -1.0, "confirm", false);
-%!
-%! Ae = [ 0.8924 0.3887 0.1285 0.1716
-%! -0.0837 0.6186 -0.6273 -0.4582
-%! 0.0052 0.1307 0.6685 -0.6755
-%! 0.0055 0.0734 -0.2148 0.4788 ];
-%!
-%! Ce = [ -0.4442 0.6663 0.3961 0.4102 ];
-%!
-%! Be = [ -0.2142
-%! -0.1968
-%! 0.0525
-%! 0.0361 ];
-%!
-%! De = [ -0.0041 ];
-%!
-%! Ke = [ -1.9513
-%! -0.1867
-%! 0.6348
-%! -0.3486 ];
-%!
-%! Qe = [ 0.0052 0.0005 -0.0017 0.0009
-%! 0.0005 0.0000 -0.0002 0.0001
-%! -0.0017 -0.0002 0.0006 -0.0003
-%! 0.0009 0.0001 -0.0003 0.0002 ];
-%!
-%! Rye = [ 0.0012 ];
-%!
-%! Se = [ -0.0025
-%! -0.0002
-%! 0.0008
-%! -0.0005 ];
-%!
-%! X0e = [ -11.496422
-%! -0.718576
-%! -0.014211
-%! 0.500073 ]; # X0e is not from SLICOT
-%!
-%! ## The SLICOT test for IB01CD uses COMUSE=C, not COMUSE=U.
-%! ## This means that they don't use the matrices B and D
-%! ## computed by IB01BD. They use only A and C from IB01BD,
-%! ## while B and D are from SLICOT routine IB01CD.
-%! ## Therefore they get slightly different matrices B and D
-%! ## and finally a different initial state vector X0.
-%!
-%!assert (SYS.A, Ae, 1e-4);
-%!assert (SYS.B, Be, 1e-4);
-%!assert (SYS.C, Ce, 1e-4);
-%!assert (SYS.D, De, 1e-4);
-%!assert (INFO.K, Ke, 1e-4);
-%!assert (INFO.Q, Qe, 1e-4);
-%!assert (INFO.Ry, Rye, 1e-4);
-%!assert (INFO.S, Se, 1e-4);
-%!assert (X0, X0e, 1e-4);
-
-
-## [96-003] Data of a 120 MW power plant (Pont-sur-Sambre, France)
-%!shared SYS, Ae, Be, Ce, De
-%! U = [ -811 -592 421 -680 -681
-%! -812 -619 477 -685 -651
-%! -817 -565 538 -678 -677
-%! -695 -725 536 -674 -702
-%! -697 -571 531 -676 -685
-%! -697 -618 533 -681 -721
-%! -702 -579 549 -677 -699
-%! -703 -487 575 -677 -694
-%! -705 -449 561 -679 -678
-%! -705 -431 563 -680 -692
-%! -707 -502 561 -679 -686
-%! -707 -583 530 -676 -751
-%! -710 -458 540 -677 -700
-%! -713 -469 543 -679 -731
-%! -715 -506 549 -684 -635
-%! -713 -590 532 -681 -715
-%! -714 -582 528 -676 -696
-%! -713 -575 538 -679 -690
-%! -716 -382 557 -682 -701
-%! -716 -476 558 -679 -690
-%! -718 -425 565 -678 -686
-%! -719 -409 562 -679 -615
-%! -719 -508 523 -677 -737
-%! -721 -569 523 -679 -722
-%! -723 -434 542 -681 -735
-%! -723 -395 544 -676 -704
-%! -723 -428 542 -677 -729
-%! -722 -402 537 -677 -706
-%! -725 -380 534 -681 -696
-%! -726 -324 549 -676 -701
-%! -726 -211 578 -675 -614
-%! -727 -113 569 -677 -738
-%! -727 -208 554 -676 -737
-%! -727 -320 523 -684 -697
-%! -727 944 605 -680 -587
-%! -729 396 654 -681 -708
-%! -729 754 637 -682 -685
-%! -727 419 522 -677 -715
-%! -729 378 494 -668 -703
-%! -728 363 493 -669 -682
-%! -729 390 496 -665 -713
-%! -729 106 483 -664 -711
-%! -729 32 495 -661 -718
-%! -729 517 585 -661 -641
-%! -729 455 625 -659 -703
-%! -730 521 649 -687 -669
-%! -730 540 627 -689 -705
-%! -731 514 605 -694 -682
-%! -585 525 558 -685 -611
-%! -586 534 520 -680 -668
-%! -586 539 531 -681 -679
-%! -585 519 507 -682 -663
-%! -588 513 505 -667 -668
-%! -587 509 539 -680 -616
-%! -587 512 535 -668 -628
-%! -588 514 557 -667 -648
-%! -588 553 563 -676 -613
-%! -589 519 559 -684 -638
-%! -589 521 563 -682 -652
-%! -588 518 547 -678 -597
-%! -589 552 549 -688 -630
-%! -589 520 535 -685 -623
-%! -589 547 542 -678 -619
-%! -589 549 531 -684 -524
-%! -588 544 522 -1540 -580
-%! -588 564 555 -1538 -584
-%! -588 684 545 -1541 -564
-%! -590 558 546 -1541 -609
-%! -589 552 537 -1550 -601
-%! -591 532 526 -1548 -580
-%! -590 544 524 -1542 -565
-%! -591 559 535 -1538 -604
-%! -592 555 542 -1548 -629
-%! -591 577 532 -1549 -587
-%! -593 581 530 -1543 -585
-%! -592 562 540 -1548 -583
-%! -591 568 546 -1536 -587
-%! -593 550 557 -1533 -569
-%! -592 550 537 -1518 -527
-%! -593 568 551 -1533 -582
-%! -590 528 540 -1529 -492
-%! -590 542 532 -1525 -585
-%! -590 556 535 -1522 -606
-%! -591 637 535 -1516 -571
-%! -591 608 539 -1512 -582
-%! -591 545 527 -1510 -577
-%! ...
[truncated message content] |
