From: <par...@us...> - 2011-10-16 09:47:10
|
Revision: 8752 http://octave.svn.sourceforge.net/octave/?rev=8752&view=rev Author: paramaniac Date: 2011-10-16 09:47:01 +0000 (Sun, 16 Oct 2011) Log Message: ----------- control: add identification routines Modified Paths: -------------- trunk/octave-forge/extra/control-devel/src/MA02AD.f trunk/octave-forge/extra/control-devel/src/MB01SD.f trunk/octave-forge/extra/control-devel/src/MB02UD.f trunk/octave-forge/extra/control-devel/src/MB03UD.f trunk/octave-forge/extra/control-devel/src/MB04OD.f trunk/octave-forge/extra/control-devel/src/MB04OY.f trunk/octave-forge/extra/control-devel/src/SB02MR.f trunk/octave-forge/extra/control-devel/src/SB02MS.f trunk/octave-forge/extra/control-devel/src/SB02MV.f trunk/octave-forge/extra/control-devel/src/SB02MW.f trunk/octave-forge/extra/control-devel/src/SB04PX.f trunk/octave-forge/extra/control-devel/src/TB01WD.f trunk/octave-forge/extra/control-devel/src/select.f Added Paths: ----------- trunk/octave-forge/extra/control-devel/DESCRIPTION trunk/octave-forge/extra/control-devel/INDEX trunk/octave-forge/extra/control-devel/devel/makefile_conred.m trunk/octave-forge/extra/control-devel/devel/makefile_ident.m trunk/octave-forge/extra/control-devel/devel/makefile_modred.m trunk/octave-forge/extra/control-devel/src/IB01AD.f trunk/octave-forge/extra/control-devel/src/IB01BD.f trunk/octave-forge/extra/control-devel/src/IB01CD.f trunk/octave-forge/extra/control-devel/src/IB01MD.f trunk/octave-forge/extra/control-devel/src/IB01MY.f trunk/octave-forge/extra/control-devel/src/IB01ND.f trunk/octave-forge/extra/control-devel/src/IB01OD.f trunk/octave-forge/extra/control-devel/src/IB01OY.f trunk/octave-forge/extra/control-devel/src/IB01PD.f trunk/octave-forge/extra/control-devel/src/IB01PX.f trunk/octave-forge/extra/control-devel/src/IB01PY.f trunk/octave-forge/extra/control-devel/src/IB01QD.f trunk/octave-forge/extra/control-devel/src/IB01RD.f trunk/octave-forge/extra/control-devel/src/MA02ED.f trunk/octave-forge/extra/control-devel/src/MA02FD.f trunk/octave-forge/extra/control-devel/src/MB01RU.f trunk/octave-forge/extra/control-devel/src/MB01RX.f trunk/octave-forge/extra/control-devel/src/MB01RY.f trunk/octave-forge/extra/control-devel/src/MB01TD.f trunk/octave-forge/extra/control-devel/src/MB01UD.f trunk/octave-forge/extra/control-devel/src/MB01VD.f trunk/octave-forge/extra/control-devel/src/MB02PD.f trunk/octave-forge/extra/control-devel/src/MB02QY.f trunk/octave-forge/extra/control-devel/src/MB03OD.f trunk/octave-forge/extra/control-devel/src/MB04ID.f trunk/octave-forge/extra/control-devel/src/MB04IY.f trunk/octave-forge/extra/control-devel/src/MB04KD.f trunk/octave-forge/extra/control-devel/src/Makefile trunk/octave-forge/extra/control-devel/src/SB02MT.f trunk/octave-forge/extra/control-devel/src/SB02ND.f trunk/octave-forge/extra/control-devel/src/SB02QD.f trunk/octave-forge/extra/control-devel/src/SB02RD.f trunk/octave-forge/extra/control-devel/src/SB02RU.f trunk/octave-forge/extra/control-devel/src/SB02SD.f trunk/octave-forge/extra/control-devel/src/SB03MV.f trunk/octave-forge/extra/control-devel/src/SB03MW.f trunk/octave-forge/extra/control-devel/src/SB03MX.f trunk/octave-forge/extra/control-devel/src/SB03MY.f trunk/octave-forge/extra/control-devel/src/SB03QX.f trunk/octave-forge/extra/control-devel/src/SB03QY.f trunk/octave-forge/extra/control-devel/src/SB03SX.f trunk/octave-forge/extra/control-devel/src/SB03SY.f Removed Paths: ------------- trunk/octave-forge/extra/control-devel/devel/makefile_slconred.m trunk/octave-forge/extra/control-devel/devel/makefile_slmodred.m Added: trunk/octave-forge/extra/control-devel/DESCRIPTION =================================================================== --- trunk/octave-forge/extra/control-devel/DESCRIPTION (rev 0) +++ trunk/octave-forge/extra/control-devel/DESCRIPTION 2011-10-16 09:47:01 UTC (rev 8752) @@ -0,0 +1,11 @@ +Name: Control-Devel +Version: 0.1.50 +Date: 2011-10-16 +Author: Lukas Reichlin <luk...@gm...> +Maintainer: Lukas Reichlin <luk...@gm...> +Title: Control Systems Developer's Playground +Description: SLICOT system identification plus model and controller reduction +Depends: octave (>= 3.3.90), control (>= 2.2.0) +Autoload: yes +License: GPL version 3 or later +Url: http://octave.sf.net, http://www.slicot.org Added: trunk/octave-forge/extra/control-devel/INDEX =================================================================== --- trunk/octave-forge/extra/control-devel/INDEX (rev 0) +++ trunk/octave-forge/extra/control-devel/INDEX 2011-10-16 09:47:01 UTC (rev 8752) @@ -0,0 +1,14 @@ +control-devel >> Control Theory +Examples +Experimental Data Handling + @iddata/get + iddata + @iddata/set + @iddata/size +System Identification + moesp + n4sid +Model Reduction + hna +Controller Reduction + Copied: trunk/octave-forge/extra/control-devel/devel/makefile_conred.m (from rev 8751, trunk/octave-forge/extra/control-devel/devel/makefile_slconred.m) =================================================================== --- trunk/octave-forge/extra/control-devel/devel/makefile_conred.m (rev 0) +++ trunk/octave-forge/extra/control-devel/devel/makefile_conred.m 2011-10-16 09:47:01 UTC (rev 8752) @@ -0,0 +1,26 @@ +homedir = pwd (); +develdir = fileparts (which ("makefile_conred")); +srcdir = [develdir, "/../src"]; +cd (srcdir); + +mkoctfile SB16AD.f TB01ID.f SB16AY.f TB01KD.f AB09IX.f \ + MB04OD.f MB01WD.f SB03OD.f MB03UD.f AB05PD.f \ + AB09DD.f AB07ND.f TB01LD.f AB05QD.f SB03OU.f \ + MA02AD.f MB03QX.f select.f MB01YD.f MB01ZD.f \ + SB03OT.f MB04OY.f MB03QD.f MB04ND.f MB03QY.f \ + SB03OR.f SB03OY.f SB04PX.f MB04NY.f SB03OV.f + +mkoctfile SB16BD.f AB09AD.f AB09BD.f SB08GD.f SB08HD.f \ + TB01ID.f AB09AX.f MA02GD.f AB09BX.f TB01WD.f \ + MA02DD.f MB03UD.f select.f AB09DD.f SB03OU.f \ + MA02AD.f SB03OT.f MB04ND.f MB04OD.f SB03OR.f \ + SB03OY.f SB04PX.f MB04NY.f MB04OY.f SB03OV.f + +mkoctfile SB16CD.f SB16CY.f AB09IX.f SB03OD.f MB02UD.f \ + AB09DD.f MA02AD.f MB03UD.f select.f SB03OU.f \ + MB01SD.f SB03OT.f MB04ND.f MB04OD.f SB03OR.f \ + SB03OY.f SB04PX.f MB04NY.f MB04OY.f SB03OV.f + +system ("rm *.o"); +cd (homedir); + Added: trunk/octave-forge/extra/control-devel/devel/makefile_ident.m =================================================================== --- trunk/octave-forge/extra/control-devel/devel/makefile_ident.m (rev 0) +++ trunk/octave-forge/extra/control-devel/devel/makefile_ident.m 2011-10-16 09:47:01 UTC (rev 8752) @@ -0,0 +1,28 @@ +homedir = pwd (); +develdir = fileparts (which ("makefile_ident")); +srcdir = [develdir, "/../src"]; +cd (srcdir); + +## preprocess the input-output data +mkoctfile IB01AD.f IB01MD.f IB01ND.f IB01OD.f IB01MY.f \ + MB04OD.f MB03UD.f MB04ID.f MA02AD.f MB03OD.f \ + MB04IY.f IB01OY.f MA02ED.f MA02FD.f MB04OY.f + +## estimating system matrices, Kalman gain, and covariances +mkoctfile IB01BD.f IB01PD.f MA02AD.f SB02MT.f SB02RD.f \ + SB02ND.f MB02UD.f MA02ED.f IB01PY.f MB03OD.f \ + MB02QY.f IB01PX.f SB02MS.f SB02RU.f SB02SD.f \ + MB01RU.f SB02QD.f SB02MV.f SB02MW.f SB02MR.f \ + MB02PD.f MB01SD.f MB04KD.f MB03UD.f MB04OD.f \ + MB04OY.f MB01VD.f select.f MB01UD.f SB03SY.f \ + MB01RX.f SB03MX.f SB03SX.f MB01RY.f SB03QY.f \ + SB03QX.f SB03MY.f SB04PX.f SB03MV.f SB03MW.f + +## estimating the initial state +mkoctfile IB01CD.f TB01WD.f IB01RD.f IB01QD.f select.f \ + MB01TD.f MA02AD.f MB04OD.f MB04OY.f MB02UD.f \ + MB03UD.f MB01SD.f + +system ("rm *.o"); +cd (homedir); + Copied: trunk/octave-forge/extra/control-devel/devel/makefile_modred.m (from rev 8751, trunk/octave-forge/extra/control-devel/devel/makefile_slmodred.m) =================================================================== --- trunk/octave-forge/extra/control-devel/devel/makefile_modred.m (rev 0) +++ trunk/octave-forge/extra/control-devel/devel/makefile_modred.m 2011-10-16 09:47:01 UTC (rev 8752) @@ -0,0 +1,35 @@ +homedir = pwd (); +develdir = fileparts (which ("makefile_modred")); +srcdir = [develdir, "/../src"]; +cd (srcdir); + +mkoctfile AB09ID.f TB01PD.f SB08DD.f TB01ID.f TB01KD.f \ + AB09IX.f AB09IY.f SB08CD.f MB04ND.f TB01XD.f \ + MB04OD.f MB01WD.f MB03UD.f AB07MD.f SB01FY.f \ + AB09DD.f TB01LD.f SB03OU.f TB01UD.f MA02AD.f \ + MA02BD.f MB03OY.f MB03QX.f MB01PD.f select.f \ + MB01YD.f MB04NY.f MB01ZD.f SB03OT.f MB04OX.f \ + MB04OY.f MB03QD.f SB03OY.f MB03QY.f MB01QD.f \ + SB03OR.f SB03OV.f SB04PX.f + +mkoctfile "-Wl,-framework" "-Wl,vecLib" \ + slab09jd.cc \ + AB09JD.f TB01ID.f TB01KD.f AB07ND.f AB09JV.f \ + AB09JW.f AB09CX.f AG07BD.f AB08MD.f AB04MD.f \ + TB01LD.f delctg.f SB04PY.f AB09AX.f AB08NX.f \ + MB01SD.f AB09JX.f MA02AD.f TB01WD.f MB03OY.f \ + MB03PY.f MA02DD.f MB03UD.f MB03QX.f select.f \ + SB04PX.f SB03OU.f MB03QD.f MB03QY.f SB03OT.f \ + MB04ND.f MB04OD.f SB03OR.f SB03OY.f MB04NY.f \ + MB04OY.f SB03OV.f + +mkoctfile AB09HD.f TB01ID.f AB04MD.f TB01KD.f AB09HY.f \ + AB09IX.f MB03UD.f SB02MD.f AB09DD.f TB01LD.f \ + SB03OU.f MA02AD.f MB03QX.f select.f SB03OT.f \ + SB02MR.f SB02MS.f MB03QD.f SB02MU.f SB02MV.f \ + SB02MW.f MB04ND.f MB04OD.f MB03QY.f SB03OR.f \ + SB03OY.f SB04PX.f MB04NY.f MB04OY.f SB03OV.f + +system ("rm *.o"); +cd (homedir); + Deleted: trunk/octave-forge/extra/control-devel/devel/makefile_slconred.m =================================================================== --- trunk/octave-forge/extra/control-devel/devel/makefile_slconred.m 2011-10-15 19:22:25 UTC (rev 8751) +++ trunk/octave-forge/extra/control-devel/devel/makefile_slconred.m 2011-10-16 09:47:01 UTC (rev 8752) @@ -1,26 +0,0 @@ -homedir = pwd (); -develdir = fileparts (which ("makefile_slconred")); -srcdir = [develdir, "/../src"]; -cd (srcdir); - -mkoctfile SB16AD.f TB01ID.f SB16AY.f TB01KD.f AB09IX.f \ - MB04OD.f MB01WD.f SB03OD.f MB03UD.f AB05PD.f \ - AB09DD.f AB07ND.f TB01LD.f AB05QD.f SB03OU.f \ - MA02AD.f MB03QX.f select.f MB01YD.f MB01ZD.f \ - SB03OT.f MB04OY.f MB03QD.f MB04ND.f MB03QY.f \ - SB03OR.f SB03OY.f SB04PX.f MB04NY.f SB03OV.f - -mkoctfile SB16BD.f AB09AD.f AB09BD.f SB08GD.f SB08HD.f \ - TB01ID.f AB09AX.f MA02GD.f AB09BX.f TB01WD.f \ - MA02DD.f MB03UD.f select.f AB09DD.f SB03OU.f \ - MA02AD.f SB03OT.f MB04ND.f MB04OD.f SB03OR.f \ - SB03OY.f SB04PX.f MB04NY.f MB04OY.f SB03OV.f - -mkoctfile SB16CD.f SB16CY.f AB09IX.f SB03OD.f MB02UD.f \ - AB09DD.f MA02AD.f MB03UD.f select.f SB03OU.f \ - MB01SD.f SB03OT.f MB04ND.f MB04OD.f SB03OR.f \ - SB03OY.f SB04PX.f MB04NY.f MB04OY.f SB03OV.f - -system ("rm *.o"); -cd (homedir); - Deleted: trunk/octave-forge/extra/control-devel/devel/makefile_slmodred.m =================================================================== --- trunk/octave-forge/extra/control-devel/devel/makefile_slmodred.m 2011-10-15 19:22:25 UTC (rev 8751) +++ trunk/octave-forge/extra/control-devel/devel/makefile_slmodred.m 2011-10-16 09:47:01 UTC (rev 8752) @@ -1,35 +0,0 @@ -homedir = pwd (); -develdir = fileparts (which ("makefile_slmodred")); -srcdir = [develdir, "/../src"]; -cd (srcdir); - -mkoctfile AB09ID.f TB01PD.f SB08DD.f TB01ID.f TB01KD.f \ - AB09IX.f AB09IY.f SB08CD.f MB04ND.f TB01XD.f \ - MB04OD.f MB01WD.f MB03UD.f AB07MD.f SB01FY.f \ - AB09DD.f TB01LD.f SB03OU.f TB01UD.f MA02AD.f \ - MA02BD.f MB03OY.f MB03QX.f MB01PD.f select.f \ - MB01YD.f MB04NY.f MB01ZD.f SB03OT.f MB04OX.f \ - MB04OY.f MB03QD.f SB03OY.f MB03QY.f MB01QD.f \ - SB03OR.f SB03OV.f SB04PX.f - -mkoctfile "-Wl,-framework" "-Wl,vecLib" \ - slab09jd.cc \ - AB09JD.f TB01ID.f TB01KD.f AB07ND.f AB09JV.f \ - AB09JW.f AB09CX.f AG07BD.f AB08MD.f AB04MD.f \ - TB01LD.f delctg.f SB04PY.f AB09AX.f AB08NX.f \ - MB01SD.f AB09JX.f MA02AD.f TB01WD.f MB03OY.f \ - MB03PY.f MA02DD.f MB03UD.f MB03QX.f select.f \ - SB04PX.f SB03OU.f MB03QD.f MB03QY.f SB03OT.f \ - MB04ND.f MB04OD.f SB03OR.f SB03OY.f MB04NY.f \ - MB04OY.f SB03OV.f - -mkoctfile AB09HD.f TB01ID.f AB04MD.f TB01KD.f AB09HY.f \ - AB09IX.f MB03UD.f SB02MD.f AB09DD.f TB01LD.f \ - SB03OU.f MA02AD.f MB03QX.f select.f SB03OT.f \ - SB02MR.f SB02MS.f MB03QD.f SB02MU.f SB02MV.f \ - SB02MW.f MB04ND.f MB04OD.f MB03QY.f SB03OR.f \ - SB03OY.f SB04PX.f MB04NY.f MB04OY.f SB03OV.f - -system ("rm *.o"); -cd (homedir); - Added: trunk/octave-forge/extra/control-devel/src/IB01AD.f =================================================================== --- trunk/octave-forge/extra/control-devel/src/IB01AD.f (rev 0) +++ trunk/octave-forge/extra/control-devel/src/IB01AD.f 2011-10-16 09:47:01 UTC (rev 8752) @@ -0,0 +1,686 @@ + SUBROUTINE IB01AD( METH, ALG, JOBD, BATCH, CONCT, CTRL, NOBR, M, + $ L, NSMP, U, LDU, Y, LDY, N, R, LDR, SV, RCOND, + $ TOL, IWORK, DWORK, LDWORK, IWARN, INFO ) +C +C SLICOT RELEASE 5.0. +C +C Copyright (c) 2002-2010 NICONET e.V. +C +C This program is free software: you can redistribute it and/or +C modify it under the terms of the GNU General Public License as +C published by the Free Software Foundation, either version 2 of +C the License, or (at your option) any later version. +C +C This program is distributed in the hope that it will be useful, +C but WITHOUT ANY WARRANTY; without even the implied warranty of +C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +C GNU General Public License for more details. +C +C You should have received a copy of the GNU General Public License +C along with this program. If not, see +C <http://www.gnu.org/licenses/>. +C +C PURPOSE +C +C To preprocess the input-output data for estimating the matrices +C of a linear time-invariant dynamical system and to find an +C estimate of the system order. The input-output data can, +C optionally, be processed sequentially. +C +C ARGUMENTS +C +C Mode Parameters +C +C METH CHARACTER*1 +C Specifies the subspace identification method to be used, +C as follows: +C = 'M': MOESP algorithm with past inputs and outputs; +C = 'N': N4SID algorithm. +C +C ALG CHARACTER*1 +C Specifies the algorithm for computing the triangular +C factor R, as follows: +C = 'C': Cholesky algorithm applied to the correlation +C matrix of the input-output data; +C = 'F': Fast QR algorithm; +C = 'Q': QR algorithm applied to the concatenated block +C Hankel matrices. +C +C JOBD CHARACTER*1 +C Specifies whether or not the matrices B and D should later +C be computed using the MOESP approach, as follows: +C = 'M': the matrices B and D should later be computed +C using the MOESP approach; +C = 'N': the matrices B and D should not be computed using +C the MOESP approach. +C This parameter is not relevant for METH = 'N'. +C +C BATCH CHARACTER*1 +C Specifies whether or not sequential data processing is to +C be used, and, for sequential processing, whether or not +C the current data block is the first block, an intermediate +C block, or the last block, as follows: +C = 'F': the first block in sequential data processing; +C = 'I': an intermediate block in sequential data +C processing; +C = 'L': the last block in sequential data processing; +C = 'O': one block only (non-sequential data processing). +C NOTE that when 100 cycles of sequential data processing +C are completed for BATCH = 'I', a warning is +C issued, to prevent for an infinite loop. +C +C CONCT CHARACTER*1 +C Specifies whether or not the successive data blocks in +C sequential data processing belong to a single experiment, +C as follows: +C = 'C': the current data block is a continuation of the +C previous data block and/or it will be continued +C by the next data block; +C = 'N': there is no connection between the current data +C block and the previous and/or the next ones. +C This parameter is not used if BATCH = 'O'. +C +C CTRL CHARACTER*1 +C Specifies whether or not the user's confirmation of the +C system order estimate is desired, as follows: +C = 'C': user's confirmation; +C = 'N': no confirmation. +C If CTRL = 'C', a reverse communication routine, IB01OY, +C is indirectly called (by SLICOT Library routine IB01OD), +C and, after inspecting the singular values and system order +C estimate, n, the user may accept n or set a new value. +C IB01OY is not called if CTRL = 'N'. +C +C Input/Output Parameters +C +C NOBR (input) INTEGER +C The number of block rows, s, in the input and output +C block Hankel matrices to be processed. NOBR > 0. +C (In the MOESP theory, NOBR should be larger than n, +C the estimated dimension of state vector.) +C +C M (input) INTEGER +C The number of system inputs. M >= 0. +C When M = 0, no system inputs are processed. +C +C L (input) INTEGER +C The number of system outputs. L > 0. +C +C NSMP (input) INTEGER +C The number of rows of matrices U and Y (number of +C samples, t). (When sequential data processing is used, +C NSMP is the number of samples of the current data +C block.) +C NSMP >= 2*(M+L+1)*NOBR - 1, for non-sequential +C processing; +C NSMP >= 2*NOBR, for sequential processing. +C The total number of samples when calling the routine with +C BATCH = 'L' should be at least 2*(M+L+1)*NOBR - 1. +C The NSMP argument may vary from a cycle to another in +C sequential data processing, but NOBR, M, and L should +C be kept constant. For efficiency, it is advisable to use +C NSMP as large as possible. +C +C U (input) DOUBLE PRECISION array, dimension (LDU,M) +C The leading NSMP-by-M part of this array must contain the +C t-by-m input-data sequence matrix U, +C U = [u_1 u_2 ... u_m]. Column j of U contains the +C NSMP values of the j-th input component for consecutive +C time increments. +C If M = 0, this array is not referenced. +C +C LDU INTEGER +C The leading dimension of the array U. +C LDU >= NSMP, if M > 0; +C LDU >= 1, if M = 0. +C +C Y (input) DOUBLE PRECISION array, dimension (LDY,L) +C The leading NSMP-by-L part of this array must contain the +C t-by-l output-data sequence matrix Y, +C Y = [y_1 y_2 ... y_l]. Column j of Y contains the +C NSMP values of the j-th output component for consecutive +C time increments. +C +C LDY INTEGER +C The leading dimension of the array Y. LDY >= NSMP. +C +C N (output) INTEGER +C The estimated order of the system. +C If CTRL = 'C', the estimated order has been reset to a +C value specified by the user. +C +C R (output or input/output) DOUBLE PRECISION array, dimension +C ( LDR,2*(M+L)*NOBR ) +C On exit, if ALG = 'C' and BATCH = 'F' or 'I', the leading +C 2*(M+L)*NOBR-by-2*(M+L)*NOBR upper triangular part of this +C array contains the current upper triangular part of the +C correlation matrix in sequential data processing. +C If ALG = 'F' and BATCH = 'F' or 'I', the array R is not +C referenced. +C On exit, if INFO = 0, ALG = 'Q', and BATCH = 'F' or 'I', +C the leading 2*(M+L)*NOBR-by-2*(M+L)*NOBR upper triangular +C part of this array contains the current upper triangular +C factor R from the QR factorization of the concatenated +C block Hankel matrices. Denote R_ij, i,j = 1:4, the +C ij submatrix of R, partitioned by M*NOBR, M*NOBR, +C L*NOBR, and L*NOBR rows and columns. +C On exit, if INFO = 0 and BATCH = 'L' or 'O', the leading +C 2*(M+L)*NOBR-by-2*(M+L)*NOBR upper triangular part of +C this array contains the matrix S, the processed upper +C triangular factor R from the QR factorization of the +C concatenated block Hankel matrices, as required by other +C subroutines. Specifically, let S_ij, i,j = 1:4, be the +C ij submatrix of S, partitioned by M*NOBR, L*NOBR, +C M*NOBR, and L*NOBR rows and columns. The submatrix +C S_22 contains the matrix of left singular vectors needed +C subsequently. Useful information is stored in S_11 and +C in the block-column S_14 : S_44. For METH = 'M' and +C JOBD = 'M', the upper triangular part of S_31 contains +C the upper triangular factor in the QR factorization of the +C matrix R_1c = [ R_12' R_22' R_11' ]', and S_12 +C contains the corresponding leading part of the transformed +C matrix R_2c = [ R_13' R_23' R_14' ]'. For METH = 'N', +C the subarray S_41 : S_43 contains the transpose of the +C matrix contained in S_14 : S_34. +C The details of the contents of R need not be known if this +C routine is followed by SLICOT Library routine IB01BD. +C On entry, if ALG = 'C', or ALG = 'Q', and BATCH = 'I' or +C 'L', the leading 2*(M+L)*NOBR-by-2*(M+L)*NOBR upper +C triangular part of this array must contain the upper +C triangular matrix R computed at the previous call of this +C routine in sequential data processing. The array R need +C not be set on entry if ALG = 'F' or if BATCH = 'F' or 'O'. +C +C LDR INTEGER +C The leading dimension of the array R. +C LDR >= MAX( 2*(M+L)*NOBR, 3*M*NOBR ), +C for METH = 'M' and JOBD = 'M'; +C LDR >= 2*(M+L)*NOBR, for METH = 'M' and JOBD = 'N' or +C for METH = 'N'. +C +C SV (output) DOUBLE PRECISION array, dimension ( L*NOBR ) +C The singular values used to estimate the system order. +C +C Tolerances +C +C RCOND DOUBLE PRECISION +C The tolerance to be used for estimating the rank of +C matrices. If the user sets RCOND > 0, the given value +C of RCOND is used as a lower bound for the reciprocal +C condition number; an m-by-n matrix whose estimated +C condition number is less than 1/RCOND is considered to +C be of full rank. If the user sets RCOND <= 0, then an +C implicitly computed, default tolerance, defined by +C RCONDEF = m*n*EPS, is used instead, where EPS is the +C relative machine precision (see LAPACK Library routine +C DLAMCH). +C This parameter is not used for METH = 'M'. +C +C TOL DOUBLE PRECISION +C Absolute tolerance used for determining an estimate of +C the system order. If TOL >= 0, the estimate is +C indicated by the index of the last singular value greater +C than or equal to TOL. (Singular values less than TOL +C are considered as zero.) When TOL = 0, an internally +C computed default value, TOL = NOBR*EPS*SV(1), is used, +C where SV(1) is the maximal singular value, and EPS is +C the relative machine precision (see LAPACK Library routine +C DLAMCH). When TOL < 0, the estimate is indicated by the +C index of the singular value that has the largest +C logarithmic gap to its successor. +C +C Workspace +C +C IWORK INTEGER array, dimension (LIWORK) +C LIWORK >= (M+L)*NOBR, if METH = 'N'; +C LIWORK >= M+L, if METH = 'M' and ALG = 'F'; +C LIWORK >= 0, if METH = 'M' and ALG = 'C' or 'Q'. +C +C DWORK DOUBLE PRECISION array, dimension (LDWORK) +C On exit, if INFO = 0, DWORK(1) returns the optimal value +C of LDWORK, and, for METH = 'N', and BATCH = 'L' or +C 'O', DWORK(2) and DWORK(3) contain the reciprocal +C condition numbers of the triangular factors of the +C matrices U_f and r_1 [6]. +C On exit, if INFO = -23, DWORK(1) returns the minimum +C value of LDWORK. +C Let +C k = 0, if CONCT = 'N' and ALG = 'C' or 'Q'; +C k = 2*NOBR-1, if CONCT = 'C' and ALG = 'C' or 'Q'; +C k = 2*NOBR*(M+L+1), if CONCT = 'N' and ALG = 'F'; +C k = 2*NOBR*(M+L+2), if CONCT = 'C' and ALG = 'F'. +C The first (M+L)*k elements of DWORK should be preserved +C during successive calls of the routine with BATCH = 'F' +C or 'I', till the final call with BATCH = 'L'. +C +C LDWORK INTEGER +C The length of the array DWORK. +C LDWORK >= (4*NOBR-2)*(M+L), if ALG = 'C', BATCH = 'F' or +C 'I' and CONCT = 'C'; +C LDWORK >= 1, if ALG = 'C', BATCH = 'F' or 'I' and +C CONCT = 'N'; +C LDWORK >= max((4*NOBR-2)*(M+L), 5*L*NOBR), if METH = 'M', +C ALG = 'C', BATCH = 'L' and CONCT = 'C'; +C LDWORK >= max((2*M-1)*NOBR, (M+L)*NOBR, 5*L*NOBR), +C if METH = 'M', JOBD = 'M', ALG = 'C', +C BATCH = 'O', or +C (BATCH = 'L' and CONCT = 'N'); +C LDWORK >= 5*L*NOBR, if METH = 'M', JOBD = 'N', ALG = 'C', +C BATCH = 'O', or +C (BATCH = 'L' and CONCT = 'N'); +C LDWORK >= 5*(M+L)*NOBR+1, if METH = 'N', ALG = 'C', and +C BATCH = 'L' or 'O'; +C LDWORK >= (M+L)*2*NOBR*(M+L+3), if ALG = 'F', +C BATCH <> 'O' and CONCT = 'C'; +C LDWORK >= (M+L)*2*NOBR*(M+L+1), if ALG = 'F', +C BATCH = 'F', 'I' and CONCT = 'N'; +C LDWORK >= (M+L)*4*NOBR*(M+L+1)+(M+L)*2*NOBR, if ALG = 'F', +C BATCH = 'L' and CONCT = 'N', or +C BATCH = 'O'; +C LDWORK >= 4*(M+L)*NOBR, if ALG = 'Q', BATCH = 'F', and +C LDR >= NS = NSMP - 2*NOBR + 1; +C LDWORK >= max(4*(M+L)*NOBR, 5*L*NOBR), if METH = 'M', +C ALG = 'Q', BATCH = 'O', and LDR >= NS; +C LDWORK >= 5*(M+L)*NOBR+1, if METH = 'N', ALG = 'Q', +C BATCH = 'O', and LDR >= NS; +C LDWORK >= 6*(M+L)*NOBR, if ALG = 'Q', (BATCH = 'F' or 'O', +C and LDR < NS), or (BATCH = 'I' or +C 'L' and CONCT = 'N'); +C LDWORK >= 4*(NOBR+1)*(M+L)*NOBR, if ALG = 'Q', BATCH = 'I' +C or 'L' and CONCT = 'C'. +C The workspace used for ALG = 'Q' is +C LDRWRK*2*(M+L)*NOBR + 4*(M+L)*NOBR, +C where LDRWRK = LDWORK/(2*(M+L)*NOBR) - 2; recommended +C value LDRWRK = NS, assuming a large enough cache size. +C For good performance, LDWORK should be larger. +C +C Warning Indicator +C +C IWARN INTEGER +C = 0: no warning; +C = 1: the number of 100 cycles in sequential data +C processing has been exhausted without signaling +C that the last block of data was get; the cycle +C counter was reinitialized; +C = 2: a fast algorithm was requested (ALG = 'C' or 'F'), +C but it failed, and the QR algorithm was then used +C (non-sequential data processing); +C = 3: all singular values were exactly zero, hence N = 0 +C (both input and output were identically zero); +C = 4: the least squares problems with coefficient matrix +C U_f, used for computing the weighted oblique +C projection (for METH = 'N'), have a rank-deficient +C coefficient matrix; +C = 5: the least squares problem with coefficient matrix +C r_1 [6], used for computing the weighted oblique +C projection (for METH = 'N'), has a rank-deficient +C coefficient matrix. +C NOTE: the values 4 and 5 of IWARN have no significance +C for the identification problem. +C +C Error Indicator +C +C INFO INTEGER +C = 0: successful exit; +C < 0: if INFO = -i, the i-th argument had an illegal +C value; +C = 1: a fast algorithm was requested (ALG = 'C', or 'F') +C in sequential data processing, but it failed; the +C routine can be repeatedly called again using the +C standard QR algorithm; +C = 2: the singular value decomposition (SVD) algorithm did +C not converge. +C +C METHOD +C +C The procedure consists in three main steps, the first step being +C performed by one of the three algorithms included. +C +C 1.a) For non-sequential data processing using QR algorithm, a +C t x 2(m+l)s matrix H is constructed, where +C +C H = [ Uf' Up' Y' ], for METH = 'M', +C s+1,2s,t 1,s,t 1,2s,t +C +C H = [ U' Y' ], for METH = 'N', +C 1,2s,t 1,2s,t +C +C and Up , Uf , U , and Y are block Hankel +C 1,s,t s+1,2s,t 1,2s,t 1,2s,t +C matrices defined in terms of the input and output data [3]. +C A QR factorization is used to compress the data. +C The fast QR algorithm uses a QR factorization which exploits +C the block-Hankel structure. Actually, the Cholesky factor of H'*H +C is computed. +C +C 1.b) For sequential data processing using QR algorithm, the QR +C decomposition is done sequentially, by updating the upper +C triangular factor R. This is also performed internally if the +C workspace is not large enough to accommodate an entire batch. +C +C 1.c) For non-sequential or sequential data processing using +C Cholesky algorithm, the correlation matrix of input-output data is +C computed (sequentially, if requested), taking advantage of the +C block Hankel structure [7]. Then, the Cholesky factor of the +C correlation matrix is found, if possible. +C +C 2) A singular value decomposition (SVD) of a certain matrix is +C then computed, which reveals the order n of the system as the +C number of "non-zero" singular values. For the MOESP approach, this +C matrix is [ R_24' R_34' ]' := R(ms+1:(2m+l)s,(2m+l)s+1:2(m+l)s), +C where R is the upper triangular factor R constructed by SLICOT +C Library routine IB01MD. For the N4SID approach, a weighted +C oblique projection is computed from the upper triangular factor R +C and its SVD is then found. +C +C 3) The singular values are compared to the given, or default TOL, +C and the estimated order n is returned, possibly after user's +C confirmation. +C +C REFERENCES +C +C [1] Verhaegen M., and Dewilde, P. +C Subspace Model Identification. Part 1: The output-error +C state-space model identification class of algorithms. +C Int. J. Control, 56, pp. 1187-1210, 1992. +C +C [2] Verhaegen M. +C Subspace Model Identification. Part 3: Analysis of the +C ordinary output-error state-space model identification +C algorithm. +C Int. J. Control, 58, pp. 555-586, 1993. +C +C [3] Verhaegen M. +C Identification of the deterministic part of MIMO state space +C models given in innovations form from input-output data. +C Automatica, Vol.30, No.1, pp.61-74, 1994. +C +C [4] Van Overschee, P., and De Moor, B. +C N4SID: Subspace Algorithms for the Identification of +C Combined Deterministic-Stochastic Systems. +C Automatica, Vol.30, No.1, pp. 75-93, 1994. +C +C [5] Peternell, K., Scherrer, W. and Deistler, M. +C Statistical Analysis of Novel Subspace Identification Methods. +C Signal Processing, 52, pp. 161-177, 1996. +C +C [6] Sima, V. +C Subspace-based Algorithms for Multivariable System +C Identification. +C Studies in Informatics and Control, 5, pp. 335-344, 1996. +C +C [7] Sima, V. +C Cholesky or QR Factorization for Data Compression in +C Subspace-based Identification ? +C Proceedings of the Second NICONET Workshop on ``Numerical +C Control Software: SLICOT, a Useful Tool in Industry'', +C December 3, 1999, INRIA Rocquencourt, France, pp. 75-80, 1999. +C +C NUMERICAL ASPECTS +C +C The implemented method is numerically stable (when QR algorithm is +C used), reliable and efficient. The fast Cholesky or QR algorithms +C are more efficient, but the accuracy could diminish by forming the +C correlation matrix. +C The most time-consuming computational step is step 1: +C 2 +C The QR algorithm needs 0(t(2(m+l)s) ) floating point operations. +C 2 3 +C The Cholesky algorithm needs 0(2t(m+l) s)+0((2(m+l)s) ) floating +C point operations. +C 2 3 2 +C The fast QR algorithm needs 0(2t(m+l) s)+0(4(m+l) s ) floating +C point operations. +C 3 +C Step 2 of the algorithm requires 0(((m+l)s) ) floating point +C operations. +C +C FURTHER COMMENTS +C +C For ALG = 'Q', BATCH = 'O' and LDR < NS, or BATCH <> 'O', the +C calculations could be rather inefficient if only minimal workspace +C (see argument LDWORK) is provided. It is advisable to provide as +C much workspace as possible. Almost optimal efficiency can be +C obtained for LDWORK = (NS+2)*(2*(M+L)*NOBR), assuming that the +C cache size is large enough to accommodate R, U, Y, and DWORK. +C +C CONTRIBUTOR +C +C V. Sima, Katholieke Universiteit Leuven, Feb. 2000. +C +C REVISIONS +C +C August 2000, March 2005. +C +C KEYWORDS +C +C Cholesky decomposition, Hankel matrix, identification methods, +C multivariable systems, QR decomposition, singular value +C decomposition. +C +C ****************************************************************** +C +C .. Scalar Arguments .. + DOUBLE PRECISION RCOND, TOL + INTEGER INFO, IWARN, L, LDR, LDU, LDWORK, LDY, M, N, + $ NOBR, NSMP + CHARACTER ALG, BATCH, CONCT, CTRL, JOBD, METH +C .. Array Arguments .. + INTEGER IWORK(*) + DOUBLE PRECISION DWORK(*), R(LDR, *), SV(*), U(LDU, *), + $ Y(LDY, *) +C .. Local Scalars .. + INTEGER IWARNL, LMNOBR, LNOBR, MAXWRK, MINWRK, MNOBR, + $ NOBR21, NR, NS, NSMPSM + LOGICAL CHALG, CONNEC, CONTRL, FIRST, FQRALG, INTERM, + $ JOBDM, LAST, MOESP, N4SID, ONEBCH, QRALG +C .. External Functions .. + LOGICAL LSAME + EXTERNAL LSAME +C .. External Subroutines .. + EXTERNAL IB01MD, IB01ND, IB01OD, XERBLA +C .. Intrinsic Functions .. + INTRINSIC MAX +C .. Save Statement .. +C MAXWRK is used to store the optimal workspace. +C NSMPSM is used to sum up the NSMP values for BATCH <> 'O'. + SAVE MAXWRK, NSMPSM +C .. +C .. Executable Statements .. +C +C Decode the scalar input parameters. +C + MOESP = LSAME( METH, 'M' ) + N4SID = LSAME( METH, 'N' ) + FQRALG = LSAME( ALG, 'F' ) + QRALG = LSAME( ALG, 'Q' ) + CHALG = LSAME( ALG, 'C' ) + JOBDM = LSAME( JOBD, 'M' ) + ONEBCH = LSAME( BATCH, 'O' ) + FIRST = LSAME( BATCH, 'F' ) .OR. ONEBCH + INTERM = LSAME( BATCH, 'I' ) + LAST = LSAME( BATCH, 'L' ) .OR. ONEBCH + CONTRL = LSAME( CTRL, 'C' ) +C + IF( .NOT.ONEBCH ) THEN + CONNEC = LSAME( CONCT, 'C' ) + ELSE + CONNEC = .FALSE. + END IF +C + MNOBR = M*NOBR + LNOBR = L*NOBR + LMNOBR = LNOBR + MNOBR + NR = LMNOBR + LMNOBR + NOBR21 = 2*NOBR - 1 + IWARN = 0 + INFO = 0 + IF( FIRST ) THEN + MAXWRK = 1 + NSMPSM = 0 + END IF + NSMPSM = NSMPSM + NSMP +C +C Check the scalar input parameters. +C + IF( .NOT.( MOESP .OR. N4SID ) ) THEN + INFO = -1 + ELSE IF( .NOT.( FQRALG .OR. QRALG .OR. CHALG ) ) THEN + INFO = -2 + ELSE IF( MOESP .AND. .NOT.( JOBDM .OR. LSAME( JOBD, 'N' ) ) ) THEN + INFO = -3 + ELSE IF( .NOT.( FIRST .OR. INTERM .OR. LAST ) ) THEN + INFO = -4 + ELSE IF( .NOT. ONEBCH ) THEN + IF( .NOT.( CONNEC .OR. LSAME( CONCT, 'N' ) ) ) + $ INFO = -5 + END IF + IF( INFO.EQ.0 ) THEN + IF( .NOT.( CONTRL .OR. LSAME( CTRL, 'N' ) ) ) THEN + INFO = -6 + ELSE IF( NOBR.LE.0 ) THEN + INFO = -7 + ELSE IF( M.LT.0 ) THEN + INFO = -8 + ELSE IF( L.LE.0 ) THEN + INFO = -9 + ELSE IF( NSMP.LT.2*NOBR .OR. + $ ( LAST .AND. NSMPSM.LT.NR+NOBR21 ) ) THEN + INFO = -10 + ELSE IF( LDU.LT.1 .OR. ( M.GT.0 .AND. LDU.LT.NSMP ) ) THEN + INFO = -12 + ELSE IF( LDY.LT.NSMP ) THEN + INFO = -14 + ELSE IF( LDR.LT.NR .OR. ( MOESP .AND. JOBDM .AND. + $ LDR.LT.3*MNOBR ) ) THEN + INFO = -17 + ELSE +C +C Compute workspace. +C (Note: Comments in the code beginning "Workspace:" describe +C the minimal amount of workspace needed at that point in the +C code, as well as the preferred amount for good performance.) +C + NS = NSMP - NOBR21 + IF ( CHALG ) THEN + IF ( .NOT.LAST ) THEN + IF ( CONNEC ) THEN + MINWRK = 2*( NR - M - L ) + ELSE + MINWRK = 1 + END IF + ELSE IF ( MOESP ) THEN + IF ( CONNEC .AND. .NOT.ONEBCH ) THEN + MINWRK = MAX( 2*( NR - M - L ), 5*LNOBR ) + ELSE + MINWRK = 5*LNOBR + IF ( JOBDM ) + $ MINWRK = MAX( 2*MNOBR - NOBR, LMNOBR, MINWRK ) + END IF + ELSE + MINWRK = 5*LMNOBR + 1 + END IF + ELSE IF ( FQRALG ) THEN + IF ( .NOT.ONEBCH .AND. CONNEC ) THEN + MINWRK = NR*( M + L + 3 ) + ELSE IF ( FIRST .OR. INTERM ) THEN + MINWRK = NR*( M + L + 1 ) + ELSE + MINWRK = 2*NR*( M + L + 1 ) + NR + END IF + ELSE + MINWRK = 2*NR + IF ( ONEBCH .AND. LDR.GE.NS ) THEN + IF ( MOESP ) THEN + MINWRK = MAX( MINWRK, 5*LNOBR ) + ELSE + MINWRK = 5*LMNOBR + 1 + END IF + END IF + IF ( FIRST ) THEN + IF ( LDR.LT.NS ) THEN + MINWRK = MINWRK + NR + END IF + ELSE + IF ( CONNEC ) THEN + MINWRK = MINWRK*( NOBR + 1 ) + ELSE + MINWRK = MINWRK + NR + END IF + END IF + END IF +C + MAXWRK = MINWRK +C + IF( LDWORK.LT.MINWRK ) THEN + INFO = -23 + DWORK( 1 ) = MINWRK + END IF + END IF + END IF +C +C Return if there are illegal arguments. +C + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'IB01AD', -INFO ) + RETURN + END IF +C +C Compress the input-output data. +C Workspace: need c*(M+L)*NOBR, where c is a constant depending +C on the algorithm and the options used +C (see SLICOT Library routine IB01MD); +C prefer larger. +C + CALL IB01MD( METH, ALG, BATCH, CONCT, NOBR, M, L, NSMP, U, LDU, Y, + $ LDY, R, LDR, IWORK, DWORK, LDWORK, IWARN, INFO ) +C + IF ( INFO.EQ.1 ) THEN +C +C Error return: A fast algorithm was requested (ALG = 'C', 'F') +C in sequential data processing, but it failed. +C + RETURN + END IF +C + MAXWRK = MAX( MAXWRK, INT( DWORK( 1 ) ) ) +C + IF ( .NOT.LAST ) THEN +C +C Return to get new data. +C + RETURN + END IF +C +C Find the singular value decomposition (SVD) giving the system +C order, and perform related preliminary calculations needed for +C computing the system matrices. +C Workspace: need max( (2*M-1)*NOBR, (M+L)*NOBR, 5*L*NOBR ), +C if METH = 'M'; +C 5*(M+L)*NOBR+1, if METH = 'N'; +C prefer larger. +C + CALL IB01ND( METH, JOBD, NOBR, M, L, R, LDR, SV, RCOND, IWORK, + $ DWORK, LDWORK, IWARNL, INFO ) + IWARN = MAX( IWARN, IWARNL ) +C + IF ( INFO.EQ.2 ) THEN +C +C Error return: the singular value decomposition (SVD) algorithm +C did not converge. +C + RETURN + END IF +C +C Estimate the system order. +C + CALL IB01OD( CTRL, NOBR, L, SV, N, TOL, IWARNL, INFO ) + IWARN = MAX( IWARN, IWARNL ) +C +C Return optimal workspace in DWORK(1). +C + DWORK( 1 ) = MAX( MAXWRK, INT( DWORK( 1 ) ) ) + RETURN +C +C *** Last line of IB01AD *** + END Property changes on: trunk/octave-forge/extra/control-devel/src/IB01AD.f ___________________________________________________________________ Added: svn:executable + * Added: trunk/octave-forge/extra/control-devel/src/IB01BD.f =================================================================== --- trunk/octave-forge/extra/control-devel/src/IB01BD.f (rev 0) +++ trunk/octave-forge/extra/control-devel/src/IB01BD.f 2011-10-16 09:47:01 UTC (rev 8752) @@ -0,0 +1,791 @@ + SUBROUTINE IB01BD( METH, JOB, JOBCK, NOBR, N, M, L, NSMPL, R, + $ LDR, A, LDA, C, LDC, B, LDB, D, LDD, Q, LDQ, + $ RY, LDRY, S, LDS, K, LDK, TOL, IWORK, DWORK, + $ LDWORK, BWORK, IWARN, INFO ) +C +C SLICOT RELEASE 5.0. +C +C Copyright (c) 2002-2010 NICONET e.V. +C +C This program is free software: you can redistribute it and/or +C modify it under the terms of the GNU General Public License as +C published by the Free Software Foundation, either version 2 of +C the License, or (at your option) any later version. +C +C This program is distributed in the hope that it will be useful, +C but WITHOUT ANY WARRANTY; without even the implied warranty of +C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +C GNU General Public License for more details. +C +C You should have received a copy of the GNU General Public License +C along with this program. If not, see +C <http://www.gnu.org/licenses/>. +C +C PURPOSE +C +C To estimate the system matrices A, C, B, and D, the noise +C covariance matrices Q, Ry, and S, and the Kalman gain matrix K +C of a linear time-invariant state space model, using the +C processed triangular factor R of the concatenated block Hankel +C matrices, provided by SLICOT Library routine IB01AD. +C +C ARGUMENTS +C +C Mode Parameters +C +C METH CHARACTER*1 +C Specifies the subspace identification method to be used, +C as follows: +C = 'M': MOESP algorithm with past inputs and outputs; +C = 'N': N4SID algorithm; +C = 'C': combined method: MOESP algorithm for finding the +C matrices A and C, and N4SID algorithm for +C finding the matrices B and D. +C +C JOB CHARACTER*1 +C Specifies which matrices should be computed, as follows: +C = 'A': compute all system matrices, A, B, C, and D; +C = 'C': compute the matrices A and C only; +C = 'B': compute the matrix B only; +C = 'D': compute the matrices B and D only. +C +C JOBCK CHARACTER*1 +C Specifies whether or not the covariance matrices and the +C Kalman gain matrix are to be computed, as follows: +C = 'C': the covariance matrices only should be computed; +C = 'K': the covariance matrices and the Kalman gain +C matrix should be computed; +C = 'N': the covariance matrices and the Kalman gain matrix +C should not be computed. +C +C Input/Output Parameters +C +C NOBR (input) INTEGER +C The number of block rows, s, in the input and output +C Hankel matrices processed by other routines. NOBR > 1. +C +C N (input) INTEGER +C The order of the system. NOBR > N > 0. +C +C M (input) INTEGER +C The number of system inputs. M >= 0. +C +C L (input) INTEGER +C The number of system outputs. L > 0. +C +C NSMPL (input) INTEGER +C If JOBCK = 'C' or 'K', the total number of samples used +C for calculating the covariance matrices. +C NSMPL >= 2*(M+L)*NOBR. +C This parameter is not meaningful if JOBCK = 'N'. +C +C R (input/workspace) DOUBLE PRECISION array, dimension +C ( LDR,2*(M+L)*NOBR ) +C On entry, the leading 2*(M+L)*NOBR-by-2*(M+L)*NOBR part +C of this array must contain the relevant data for the MOESP +C or N4SID algorithms, as constructed by SLICOT Library +C routine IB01AD. Let R_ij, i,j = 1:4, be the +C ij submatrix of R (denoted S in IB01AD), partitioned +C by M*NOBR, L*NOBR, M*NOBR, and L*NOBR rows and +C columns. The submatrix R_22 contains the matrix of left +C singular vectors used. Also needed, for METH = 'N' or +C JOBCK <> 'N', are the submatrices R_11, R_14 : R_44, +C and, for METH = 'M' or 'C' and JOB <> 'C', the +C submatrices R_31 and R_12, containing the processed +C matrices R_1c and R_2c, respectively, as returned by +C SLICOT Library routine IB01AD. +C Moreover, if METH = 'N' and JOB = 'A' or 'C', the +C block-row R_41 : R_43 must contain the transpose of the +C block-column R_14 : R_34 as returned by SLICOT Library +C routine IB01AD. +C The remaining part of R is used as workspace. +C On exit, part of this array is overwritten. Specifically, +C if METH = 'M', R_22 and R_31 are overwritten if +C JOB = 'B' or 'D', and R_12, R_22, R_14 : R_34, +C and possibly R_11 are overwritten if JOBCK <> 'N'; +C if METH = 'N', all needed submatrices are overwritten. +C The details of the contents of R need not be known if +C this routine is called once just after calling the SLICOT +C Library routine IB01AD. +C +C LDR INTEGER +C The leading dimension of the array R. +C LDR >= 2*(M+L)*NOBR. +C +C A (input or output) DOUBLE PRECISION array, dimension +C (LDA,N) +C On entry, if METH = 'N' or 'C' and JOB = 'B' or 'D', +C the leading N-by-N part of this array must contain the +C system state matrix. +C If METH = 'M' or (METH = 'N' or 'C' and JOB = 'A' +C or 'C'), this array need not be set on input. +C On exit, if JOB = 'A' or 'C' and INFO = 0, the +C leading N-by-N part of this array contains the system +C state matrix. +C +C LDA INTEGER +C The leading dimension of the array A. +C LDA >= N, if JOB = 'A' or 'C', or METH = 'N' or 'C' +C and JOB = 'B' or 'D'; +C LDA >= 1, otherwise. +C +C C (input or output) DOUBLE PRECISION array, dimension +C (LDC,N) +C On entry, if METH = 'N' or 'C' and JOB = 'B' or 'D', +C the leading L-by-N part of this array must contain the +C system output matrix. +C If METH = 'M' or (METH = 'N' or 'C' and JOB = 'A' +C or 'C'), this array need not be set on input. +C On exit, if JOB = 'A' or 'C' and INFO = 0, or +C INFO = 3 (or INFO >= 0, for METH = 'M'), the leading +C L-by-N part of this array contains the system output +C matrix. +C +C LDC INTEGER +C The leading dimension of the array C. +C LDC >= L, if JOB = 'A' or 'C', or METH = 'N' or 'C' +C and JOB = 'B' or 'D'; +C LDC >= 1, otherwise. +C +C B (output) DOUBLE PRECISION array, dimension (LDB,M) +C If M > 0, JOB = 'A', 'B', or 'D' and INFO = 0, the +C leading N-by-M part of this array contains the system +C input matrix. If M = 0 or JOB = 'C', this array is +C not referenced. +C +C LDB INTEGER +C The leading dimension of the array B. +C LDB >= N, if M > 0 and JOB = 'A', 'B', or 'D'; +C LDB >= 1, if M = 0 or JOB = 'C'. +C +C D (output) DOUBLE PRECISION array, dimension (LDD,M) +C If M > 0, JOB = 'A' or 'D' and INFO = 0, the leading +C L-by-M part of this array contains the system input-output +C matrix. If M = 0 or JOB = 'C' or 'B', this array is +C not referenced. +C +C LDD INTEGER +C The leading dimension of the array D. +C LDD >= L, if M > 0 and JOB = 'A' or 'D'; +C LDD >= 1, if M = 0 or JOB = 'C' or 'B'. +C +C Q (output) DOUBLE PRECISION array, dimension (LDQ,N) +C If JOBCK = 'C' or 'K', the leading N-by-N part of this +C array contains the positive semidefinite state covariance +C matrix. If JOBCK = 'K', this matrix has been used as +C state weighting matrix for computing the Kalman gain. +C This parameter is not referenced if JOBCK = 'N'. +C +C LDQ INTEGER +C The leading dimension of the array Q. +C LDQ >= N, if JOBCK = 'C' or 'K'; +C LDQ >= 1, if JOBCK = 'N'. +C +C RY (output) DOUBLE PRECISION array, dimension (LDRY,L) +C If JOBCK = 'C' or 'K', the leading L-by-L part of this +C array contains the positive (semi)definite output +C covariance matrix. If JOBCK = 'K', this matrix has been +C used as output weighting matrix for computing the Kalman +C gain. +C This parameter is not referenced if JOBCK = 'N'. +C +C LDRY INTEGER +C The leading dimension of the array RY. +C LDRY >= L, if JOBCK = 'C' or 'K'; +C LDRY >= 1, if JOBCK = 'N'. +C +C S (output) DOUBLE PRECISION array, dimension (LDS,L) +C If JOBCK = 'C' or 'K', the leading N-by-L part of this +C array contains the state-output cross-covariance matrix. +C If JOBCK = 'K', this matrix has been used as state- +C output weighting matrix for computing the Kalman gain. +C This parameter is not referenced if JOBCK = 'N'. +C +C LDS INTEGER +C The leading dimension of the array S. +C LDS >= N, if JOBCK = 'C' or 'K'; +C LDS >= 1, if JOBCK = 'N'. +C +C K (output) DOUBLE PRECISION array, dimension ( LDK,L ) +C If JOBCK = 'K', the leading N-by-L part of this array +C contains the estimated Kalman gain matrix. +C If JOBCK = 'C' or 'N', this array is not referenced. +C +C LDK INTEGER +C The leading dimension of the array K. +C LDK >= N, if JOBCK = 'K'; +C LDK >= 1, if JOBCK = 'C' or 'N'. +C +C Tolerances +C +C TOL DOUBLE PRECISION +C The tolerance to be used for estimating the rank of +C matrices. If the user sets TOL > 0, then the given value +C of TOL is used as a lower bound for the reciprocal +C condition number; an m-by-n matrix whose estimated +C condition number is less than 1/TOL is considered to +C be of full rank. If the user sets TOL <= 0, then an +C implicitly computed, default tolerance, defined by +C TOLDEF = m*n*EPS, is used instead, where EPS is the +C relative machine precision (see LAPACK Library routine +C DLAMCH). +C +C Workspace +C +C IWORK INTEGER array, dimension (LIWORK) +C LIWORK >= max(LIW1,LIW2), where +C LIW1 = N, if METH <> 'N' and M = 0 +C or JOB = 'C' and JOBCK = 'N'; +C LIW1 = M*NOBR+N, if METH <> 'N', JOB = 'C', +C and JOBCK <> 'N'; +C LIW1 = max(L*NOBR,M*NOBR), if METH = 'M', JOB <> 'C', +C and JOBCK = 'N'; +C LIW1 = max(L*NOBR,M*NOBR+N), if METH = 'M', JOB <> 'C', +C and JOBCK = 'C' or 'K'; +C LIW1 = max(M*NOBR+N,M*(N+L)), if METH = 'N', or METH = 'C' +C and JOB <> 'C'; +C LIW2 = 0, if JOBCK <> 'K'; +C LIW2 = N*N, if JOBCK = 'K'. +C +C DWORK DOUBLE PRECISION array, dimension (LDWORK) +C On exit, if INFO = 0, DWORK(1) returns the optimal value +C of LDWORK, and DWORK(2), DWORK(3), DWORK(4), and +C DWORK(5) contain the reciprocal condition numbers of the +C triangular factors of the following matrices (defined in +C SLICOT Library routine IB01PD and in the lower level +C routines): +C GaL (GaL = Un(1:(s-1)*L,1:n)), +C R_1c (if METH = 'M' or 'C'), +C M (if JOBCK = 'C' or 'K' or METH = 'N'), and +C Q or T (see SLICOT Library routine IB01PY or IB01PX), +C respectively. +C If METH = 'N', DWORK(3) is set to one without any +C calculations. Similarly, if METH = 'M' and JOBCK = 'N', +C DWORK(4) is set to one. If M = 0 or JOB = 'C', +C DWORK(3) and DWORK(5) are set to one. +C If JOBCK = 'K' and INFO = 0, DWORK(6) to DWORK(13) +C contain information about the accuracy of the results when +C computing the Kalman gain matrix, as follows: +C DWORK(6) - reciprocal condition number of the matrix +C U11 of the Nth order system of algebraic +C equations from which the solution matrix X +C of the Riccati equation is obtained; +C DWORK(7) - reciprocal pivot growth factor for the LU +C factorization of the matrix U11; +C DWORK(8) - reciprocal condition number of the matrix +C As = A - S*inv(Ry)*C, which is inverted by +C the standard Riccati solver; +C DWORK(9) - reciprocal pivot growth factor for the LU +C factorization of the matrix As; +C DWORK(10) - reciprocal condition number of the matrix +C Ry; +C DWORK(11) - reciprocal condition number of the matrix +C Ry + C*X*C'; +C DWORK(12) - reciprocal condition number for the Riccati +C equation solution; +C DWORK(13) - forward error bound for the Riccati +C equation solution. +C On exit, if INFO = -30, DWORK(1) returns the minimum +C value of LDWORK. +C +C LDWORK INTEGER +C The length of the array DWORK. +C LDWORK >= max( LDW1,LDW2,LDW3 ), where, if METH = 'M', +C LDW1 >= max( 2*(L*NOBR-L)*N+2*N, (L*NOBR-L)*N+N*N+7*N ), +C if JOB = 'C' or JOB = 'A' and M = 0; +C LDW1 >= max( 2*(L*NOBR-L)*N+N*N+7*N, +C (L*NOBR-L)*N+N+6*M*NOBR, (L*NOBR-L)*N+N+ +C max( L+M*NOBR, L*NOBR + +C max( 3*L*NOBR+1, M ) ) ), +C if M > 0 and JOB = 'A', 'B', or 'D'; +C LDW2 >= 0, if JOBCK = 'N'; +C LDW2 >= L*NOBR*N+ +C max( (L*NOBR-L)*N+Aw+2*N+max(5*N,(2*M+L)*NOBR+L), +C 4*(M*NOBR+N)+1, M*NOBR+2*N+L ), +C if JOBCK = 'C' or 'K', +C where Aw = N+N*N, if M = 0 or JOB = 'C'; +C Aw = 0, otherwise; +C if METH = 'N', +C LDW1 >= L*NOBR*N+max( (L*NOBR-L)*N+2*N+(2*M+L)*NOBR+L, +C 2*(L*NOBR-L)*N+N*N+8*N, +C N+4*(M*NOBR+N)+1, M*NOBR+3*N+L ); +C LDW2 >= 0, if M = 0 or JOB = 'C'; +C LDW2 >= L*NOBR*N+M*NOBR*(N+L)*(M*(N+L)+1)+ +C max( (N+L)**2, 4*M*(N+L)+1 ), +C if M > 0 and JOB = 'A', 'B', or 'D'; +C and, if METH = 'C', LDW1 as +C max( LDW1 for METH = 'M', JOB = 'C', LDW1 for METH = 'N'), +C and LDW2 for METH = 'N' are used; +C LDW3 >= 0, if JOBCK <> 'K'; +C LDW3 >= max( 4*N*N+2*N*L+L*L+max( 3*L,N*L ), +C 14*N*N+12*N+5 ), if JOBCK = 'K'. +C For good performance, LDWORK should be larger. +C +C BWORK LOGICAL array, dimension ... [truncated message content] |
From: <par...@us...> - 2011-10-19 16:38:42
|
Revision: 8789 http://octave.svn.sourceforge.net/octave/?rev=8789&view=rev Author: paramaniac Date: 2011-10-19 16:38:34 +0000 (Wed, 19 Oct 2011) Log Message: ----------- control-devel: add code for fitfrd Modified Paths: -------------- trunk/octave-forge/extra/control-devel/devel/makefile_ident.m Added Paths: ----------- trunk/octave-forge/extra/control-devel/src/DG01MD.f trunk/octave-forge/extra/control-devel/src/MC01PD.f trunk/octave-forge/extra/control-devel/src/SB10YD.f trunk/octave-forge/extra/control-devel/src/SB10ZP.f trunk/octave-forge/extra/control-devel/src/TD03AY.f trunk/octave-forge/extra/control-devel/src/TD04AD.f Modified: trunk/octave-forge/extra/control-devel/devel/makefile_ident.m =================================================================== --- trunk/octave-forge/extra/control-devel/devel/makefile_ident.m 2011-10-19 15:05:10 UTC (rev 8788) +++ trunk/octave-forge/extra/control-devel/devel/makefile_ident.m 2011-10-19 16:38:34 UTC (rev 8789) @@ -23,6 +23,12 @@ MB01TD.f MA02AD.f MB04OD.f MB04OY.f MB02UD.f \ MB03UD.f MB01SD.f +## fit state-space model to frequency response data +mkoctfile SB10YD.f DG01MD.f AB04MD.f SB10ZP.f AB07ND.f \ + MC01PD.f TD04AD.f TD03AY.f TB01PD.f TB01XD.f \ + AB07MD.f TB01UD.f TB01ID.f MB01PD.f MB03OY.f \ + MB01QD.f + system ("rm *.o"); cd (homedir); Added: trunk/octave-forge/extra/control-devel/src/DG01MD.f =================================================================== --- trunk/octave-forge/extra/control-devel/src/DG01MD.f (rev 0) +++ trunk/octave-forge/extra/control-devel/src/DG01MD.f 2011-10-19 16:38:34 UTC (rev 8789) @@ -0,0 +1,235 @@ + SUBROUTINE DG01MD( INDI, N, XR, XI, INFO ) +C +C SLICOT RELEASE 5.0. +C +C Copyright (c) 2002-2010 NICONET e.V. +C +C This program is free software: you can redistribute it and/or +C modify it under the terms of the GNU General Public License as +C published by the Free Software Foundation, either version 2 of +C the License, or (at your option) any later version. +C +C This program is distributed in the hope that it will be useful, +C but WITHOUT ANY WARRANTY; without even the implied warranty of +C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +C GNU General Public License for more details. +C +C You should have received a copy of the GNU General Public License +C along with this program. If not, see +C <http://www.gnu.org/licenses/>. +C +C PURPOSE +C +C To compute the discrete Fourier transform, or inverse transform, +C of a complex signal. +C +C ARGUMENTS +C +C Mode Parameters +C +C INDI CHARACTER*1 +C Indicates whether a Fourier transform or inverse Fourier +C transform is to be performed as follows: +C = 'D': (Direct) Fourier transform; +C = 'I': Inverse Fourier transform. +C +C Input/Output Parameters +C +C N (input) INTEGER +C The number of complex samples. N must be a power of 2. +C N >= 2. +C +C XR (input/output) DOUBLE PRECISION array, dimension (N) +C On entry, this array must contain the real part of either +C the complex signal z if INDI = 'D', or f(z) if INDI = 'I'. +C On exit, this array contains either the real part of the +C computed Fourier transform f(z) if INDI = 'D', or the +C inverse Fourier transform z of f(z) if INDI = 'I'. +C +C XI (input/output) DOUBLE PRECISION array, dimension (N) +C On entry, this array must contain the imaginary part of +C either z if INDI = 'D', or f(z) if INDI = 'I'. +C On exit, this array contains either the imaginary part of +C f(z) if INDI = 'D', or z if INDI = 'I'. +C +C Error Indicator +C +C INFO INTEGER +C = 0: successful exit; +C < 0: if INFO = -i, the i-th argument had an illegal +C value. +C +C METHOD +C +C If INDI = 'D', then the routine performs a discrete Fourier +C transform on the complex signal Z(i), i = 1,2,...,N. If the result +C is denoted by FZ(k), k = 1,2,...,N, then the relationship between +C Z and FZ is given by the formula: +C +C N ((k-1)*(i-1)) +C FZ(k) = SUM ( Z(i) * V ), +C i=1 +C 2 +C where V = exp( -2*pi*j/N ) and j = -1. +C +C If INDI = 'I', then the routine performs an inverse discrete +C Fourier transform on the complex signal FZ(k), k = 1,2,...,N. If +C the result is denoted by Z(i), i = 1,2,...,N, then the +C relationship between Z and FZ is given by the formula: +C +C N ((k-1)*(i-1)) +C Z(i) = SUM ( FZ(k) * W ), +C k=1 +C +C where W = exp( 2*pi*j/N ). +C +C Note that a discrete Fourier transform, followed by an inverse +C discrete Fourier transform, will result in a signal which is a +C factor N larger than the original input signal. +C +C REFERENCES +C +C [1] Rabiner, L.R. and Rader, C.M. +C Digital Signal Processing. +C IEEE Press, 1972. +C +C NUMERICAL ASPECTS +C +C The algorithm requires 0( N*log(N) ) operations. +C +C CONTRIBUTOR +C +C Release 3.0: V. Sima, Katholieke Univ. Leuven, Belgium, Feb. 1997. +C Supersedes Release 2.0 routine DG01AD by R. Dekeyser, State +C University of Gent, Belgium. +C +C REVISIONS +C +C - +C +C KEYWORDS +C +C Complex signals, digital signal processing, fast Fourier +C transform. +C +C ****************************************************************** +C +C .. Parameters .. + DOUBLE PRECISION ZERO, HALF, ONE, TWO, EIGHT + PARAMETER ( ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0, + $ TWO = 2.0D0, EIGHT = 8.0D0 ) +C .. Scalar Arguments .. + CHARACTER INDI + INTEGER INFO, N +C .. Array Arguments .. + DOUBLE PRECISION XI(*), XR(*) +C .. Local Scalars .. + LOGICAL LINDI + INTEGER I, J, K, L, M + DOUBLE PRECISION PI2, TI, TR, WHELP, WI, WR, WSTPI, WSTPR +C .. External Functions .. + LOGICAL LSAME + EXTERNAL LSAME +C .. External Subroutines .. + EXTERNAL XERBLA +C .. Intrinsic Functions .. + INTRINSIC ATAN, DBLE, MOD, SIN +C .. Executable Statements .. +C + INFO = 0 + LINDI = LSAME( INDI, 'D' ) +C +C Test the input scalar arguments. +C + IF( .NOT.LINDI .AND. .NOT.LSAME( INDI, 'I' ) ) THEN + INFO = -1 + ELSE + J = 0 + IF( N.GE.2 ) THEN + J = N +C WHILE ( MOD( J, 2 ).EQ.0 ) DO + 10 CONTINUE + IF ( MOD( J, 2 ).EQ.0 ) THEN + J = J/2 + GO TO 10 + END IF +C END WHILE 10 + END IF + IF ( J.NE.1 ) INFO = -2 + END IF +C + IF ( INFO.NE.0 ) THEN +C +C Error return. +C + CALL XERBLA( 'DG01MD', -INFO ) + RETURN + END IF +C +C Inplace shuffling of data. +C + J = 1 +C + DO 30 I = 1, N + IF ( J.GT.I ) THEN + TR = XR(I) + TI = XI(I) + XR(I) = XR(J) + XI(I) = XI(J) + XR(J) = TR + XI(J) = TI + END IF + K = N/2 +C REPEAT + 20 IF ( J.GT.K ) THEN + J = J - K + K = K/2 + IF ( K.GE.2 ) GO TO 20 + END IF +C UNTIL ( K.LT.2 ) + J = J + K + 30 CONTINUE +C +C Transform by decimation in time. +C + PI2 = EIGHT*ATAN( ONE ) + IF ( LINDI ) PI2 = -PI2 +C + I = 1 +C +C WHILE ( I.LT.N ) DO +C + 40 IF ( I.LT.N ) THEN + L = 2*I + WHELP = PI2/DBLE( L ) + WSTPI = SIN( WHELP ) + WHELP = SIN( HALF*WHELP ) + WSTPR = -TWO*WHELP*WHELP + WR = ONE + WI = ZERO +C + DO 60 J = 1, I +C + DO 50 K = J, N, L + M = K + I + TR = WR*XR(M) - WI*XI(M) + TI = WR*XI(M) + WI*XR(M) + XR(M) = XR(K) - TR + XI(M) = XI(K) - TI + XR(K) = XR(K) + TR + XI(K) = XI(K) + TI + 50 CONTINUE +C + WHELP = WR + WR = WR + WR*WSTPR - WI*WSTPI + WI = WI + WHELP*WSTPI + WI*WSTPR + 60 CONTINUE +C + I = L + GO TO 40 +C END WHILE 40 + END IF +C + RETURN +C *** Last line of DG01MD *** + END Property changes on: trunk/octave-forge/extra/control-devel/src/DG01MD.f ___________________________________________________________________ Added: svn:executable + * Added: trunk/octave-forge/extra/control-devel/src/MC01PD.f =================================================================== --- trunk/octave-forge/extra/control-devel/src/MC01PD.f (rev 0) +++ trunk/octave-forge/extra/control-devel/src/MC01PD.f 2011-10-19 16:38:34 UTC (rev 8789) @@ -0,0 +1,159 @@ + SUBROUTINE MC01PD( K, REZ, IMZ, P, DWORK, INFO ) +C +C SLICOT RELEASE 5.0. +C +C Copyright (c) 2002-2010 NICONET e.V. +C +C This program is free software: you can redistribute it and/or +C modify it under the terms of the GNU General Public License as +C published by the Free Software Foundation, either version 2 of +C the License, or (at your option) any later version. +C +C This program is distributed in the hope that it will be useful, +C but WITHOUT ANY WARRANTY; without even the implied warranty of +C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +C GNU General Public License for more details. +C +C You should have received a copy of the GNU General Public License +C along with this program. If not, see +C <http://www.gnu.org/licenses/>. +C +C PURPOSE +C +C To compute the coefficients of a real polynomial P(x) from its +C zeros. +C +C ARGUMENTS +C +C Input/Output Parameters +C +C K (input) INTEGER +C The number of zeros (and hence the degree) of P(x). +C K >= 0. +C +C REZ (input) DOUBLE PRECISION array, dimension (K) +C IMZ (input) DOUBLE PRECISION array, dimension (K) +C The real and imaginary parts of the i-th zero of P(x) +C must be stored in REZ(i) and IMZ(i), respectively, where +C i = 1, 2, ..., K. The zeros may be supplied in any order, +C except that complex conjugate zeros must appear +C consecutively. +C +C P (output) DOUBLE PRECISION array, dimension (K+1) +C This array contains the coefficients of P(x) in increasing +C powers of x. If K = 0, then P(1) is set to one. +C +C Workspace +C +C DWORK DOUBLE PRECISION array, dimension (K+1) +C If K = 0, this array is not referenced. +C +C Error Indicator +C +C INFO INTEGER +C = 0: successful exit; +C < 0: if INFO = -i, the i-th argument had an illegal +C value; +C > 0: if INFO = i, (REZ(i),IMZ(i)) is a complex zero but +C (REZ(i-1),IMZ(i-1)) is not its conjugate. +C +C METHOD +C +C The routine computes the coefficients of the real K-th degree +C polynomial P(x) as +C +C P(x) = (x - r(1)) * (x - r(2)) * ... * (x - r(K)) +C +C where r(i) = (REZ(i),IMZ(i)). +C +C Note that REZ(i) = REZ(j) and IMZ(i) = -IMZ(j) if r(i) and r(j) +C form a complex conjugate pair (where i <> j), and that IMZ(i) = 0 +C if r(i) is real. +C +C NUMERICAL ASPECTS +C +C None. +C +C CONTRIBUTOR +C +C Release 3.0: V. Sima, Katholieke Univ. Leuven, Belgium, Mar. 1997. +C Supersedes Release 2.0 routine MC01DD by A.J. Geurts. +C +C REVISIONS +C +C V. Sima, May 2002. +C +C KEYWORDS +C +C Elementary polynomial operations, polynomial operations. +C +C ****************************************************************** +C +C .. Parameters .. + DOUBLE PRECISION ZERO, ONE + PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) +C .. Scalar Arguments .. + INTEGER INFO, K +C .. Array Arguments .. + DOUBLE PRECISION DWORK(*), IMZ(*), P(*), REZ(*) +C .. Local Scalars .. + INTEGER I + DOUBLE PRECISION U, V +C .. External Subroutines .. + EXTERNAL DAXPY, DCOPY, XERBLA +C .. Executable Statements .. +C +C Test the input scalar arguments. +C + IF( K.LT.0 ) THEN + INFO = -1 +C +C Error return. +C + CALL XERBLA( 'MC01PD', -INFO ) + RETURN + END IF +C +C Quick return if possible. +C + INFO = 0 + P(1) = ONE + IF ( K.EQ.0 ) + $ RETURN +C + I = 1 +C WHILE ( I <= K ) DO + 20 IF ( I.LE.K ) THEN + U = REZ(I) + V = IMZ(I) + DWORK(1) = ZERO +C + IF ( V.EQ.ZERO ) THEN + CALL DCOPY( I, P, 1, DWORK(2), 1 ) + CALL DAXPY( I, -U, P, 1, DWORK, 1 ) + I = I + 1 +C + ELSE + IF ( I.EQ.K ) THEN + INFO = K + RETURN + ELSE IF ( ( U.NE.REZ(I+1) ) .OR. ( V.NE.-IMZ(I+1) ) ) THEN + INFO = I + 1 + RETURN + END IF +C + DWORK(2) = ZERO + CALL DCOPY( I, P, 1, DWORK(3), 1 ) + CALL DAXPY( I, -(U + U), P, 1, DWORK(2), 1 ) + CALL DAXPY( I, U**2+V**2, P, 1, DWORK, 1 ) + I = I + 2 + END IF +C + CALL DCOPY( I, DWORK, 1, P, 1 ) + GO TO 20 + END IF +C END WHILE 20 +C + RETURN +C *** Last line of MC01PD *** + END Property changes on: trunk/octave-forge/extra/control-devel/src/MC01PD.f ___________________________________________________________________ Added: svn:executable + * Added: trunk/octave-forge/extra/control-devel/src/SB10YD.f =================================================================== --- trunk/octave-forge/extra/control-devel/src/SB10YD.f (rev 0) +++ trunk/octave-forge/extra/control-devel/src/SB10YD.f 2011-10-19 16:38:34 UTC (rev 8789) @@ -0,0 +1,689 @@ + SUBROUTINE SB10YD( DISCFL, FLAG, LENDAT, RFRDAT, IFRDAT, OMEGA, N, + $ A, LDA, B, C, D, TOL, IWORK, DWORK, LDWORK, + $ ZWORK, LZWORK, INFO ) +C +C SLICOT RELEASE 5.0. +C +C Copyright (c) 2002-2010 NICONET e.V. +C +C This program is free software: you can redistribute it and/or +C modify it under the terms of the GNU General Public License as +C published by the Free Software Foundation, either version 2 of +C the License, or (at your option) any later version. +C +C This program is distributed in the hope that it will be useful, +C but WITHOUT ANY WARRANTY; without even the implied warranty of +C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +C GNU General Public License for more details. +C +C You should have received a copy of the GNU General Public License +C along with this program. If not, see +C <http://www.gnu.org/licenses/>. +C +C PURPOSE +C +C To fit a supplied frequency response data with a stable, minimum +C phase SISO (single-input single-output) system represented by its +C matrices A, B, C, D. It handles both discrete- and continuous-time +C cases. +C +C ARGUMENTS +C +C Input/Output parameters +C +C DISCFL (input) INTEGER +C Indicates the type of the system, as follows: +C = 0: continuous-time system; +C = 1: discrete-time system. +C +C FLAG (input) INTEGER +C If FLAG = 0, then the system zeros and poles are not +C constrained. +C If FLAG = 1, then the system zeros and poles will have +C negative real parts in the continuous-time case, or moduli +C less than 1 in the discrete-time case. Consequently, FLAG +C must be equal to 1 in mu-synthesis routines. +C +C LENDAT (input) INTEGER +C The length of the vectors RFRDAT, IFRDAT and OMEGA. +C LENDAT >= 2. +C +C RFRDAT (input) DOUBLE PRECISION array, dimension (LENDAT) +C The real part of the frequency data to be fitted. +C +C IFRDAT (input) DOUBLE PRECISION array, dimension (LENDAT) +C The imaginary part of the frequency data to be fitted. +C +C OMEGA (input) DOUBLE PRECISION array, dimension (LENDAT) +C The frequencies corresponding to RFRDAT and IFRDAT. +C These values must be nonnegative and monotonically +C increasing. Additionally, for discrete-time systems +C they must be between 0 and PI. +C +C N (input/output) INTEGER +C On entry, the desired order of the system to be fitted. +C N <= LENDAT-1. +C On exit, the order of the obtained system. The value of N +C could only be modified if N > 0 and FLAG = 1. +C +C A (output) DOUBLE PRECISION array, dimension (LDA,N) +C The leading N-by-N part of this array contains the +C matrix A. If FLAG = 1, then A is in an upper Hessenberg +C form, and corresponds to a minimal realization. +C +C LDA INTEGER +C The leading dimension of the array A. LDA >= MAX(1,N). +C +C B (output) DOUBLE PRECISION array, dimension (N) +C The computed vector B. +C +C C (output) DOUBLE PRECISION array, dimension (N) +C The computed vector C. If FLAG = 1, the first N-1 elements +C are zero (for the exit value of N). +C +C D (output) DOUBLE PRECISION array, dimension (1) +C The computed scalar D. +C +C Tolerances +C +C TOL DOUBLE PRECISION +C The tolerance to be used for determining the effective +C rank of matrices. If the user sets TOL > 0, then the given +C value of TOL is used as a lower bound for the reciprocal +C condition number; a (sub)matrix whose estimated condition +C number is less than 1/TOL is considered to be of full +C rank. If the user sets TOL <= 0, then an implicitly +C computed, default tolerance, defined by TOLDEF = SIZE*EPS, +C is used instead, where SIZE is the product of the matrix +C dimensions, and EPS is the machine precision (see LAPACK +C Library routine DLAMCH). +C +C Workspace +C +C IWORK INTEGER array, dimension max(2,2*N+1) +C +C DWORK DOUBLE PRECISION array, dimension (LDWORK) +C On exit, if INFO = 0, DWORK(1) returns the optimal value +C of LDWORK and DWORK(2) contains the optimal value of +C LZWORK. +C +C LDWORK INTEGER +C The length of the array DWORK. +C LDWORK = max( 2, LW1, LW2, LW3, LW4 ), where +C LW1 = 2*LENDAT + 4*HNPTS; HNPTS = 2048; +C LW2 = LENDAT + 6*HNPTS; +C MN = min( 2*LENDAT, 2*N+1 ) +C LW3 = 2*LENDAT*(2*N+1) + max( 2*LENDAT, 2*N+1 ) + +C max( MN + 6*N + 4, 2*MN + 1 ), if N > 0; +C LW3 = 4*LENDAT + 5 , if N = 0; +C LW4 = max( N*N + 5*N, 6*N + 1 + min( 1,N ) ), if FLAG = 1; +C LW4 = 0, if FLAG = 0. +C For optimum performance LDWORK should be larger. +C +C ZWORK COMPLEX*16 array, dimension (LZWORK) +C +C LZWORK INTEGER +C The length of the array ZWORK. +C LZWORK = LENDAT*(2*N+3), if N > 0; +C LZWORK = LENDAT, if N = 0. +C +C Error indicator +C +C INFO INTEGER +C = 0: successful exit; +C < 0: if INFO = -i, the i-th argument had an illegal +C value; +C = 1: if the discrete --> continuous transformation cannot +C be made; +C = 2: if the system poles cannot be found; +C = 3: if the inverse system cannot be found, i.e., D is +C (close to) zero; +C = 4: if the system zeros cannot be found; +C = 5: if the state-space representation of the new +C transfer function T(s) cannot be found; +C = 6: if the continuous --> discrete transformation cannot +C be made. +C +C METHOD +C +C First, if the given frequency data are corresponding to a +C continuous-time system, they are changed to a discrete-time +C system using a bilinear transformation with a scaled alpha. +C Then, the magnitude is obtained from the supplied data. +C Then, the frequency data are linearly interpolated around +C the unit-disc. +C Then, Oppenheim and Schafer complex cepstrum method is applied +C to get frequency data corresponding to a stable, minimum- +C phase system. This is done in the following steps: +C - Obtain LOG (magnitude) +C - Obtain IFFT of the result (DG01MD SLICOT subroutine); +C - halve the data at 0; +C - Obtain FFT of the halved data (DG01MD SLICOT subroutine); +C - Obtain EXP of the result. +C Then, the new frequency data are interpolated back to the +C original frequency. +C Then, based on these newly obtained data, the system matrices +C A, B, C, D are constructed; the very identification is +C performed by Least Squares Method using DGELSY LAPACK subroutine. +C If needed, a discrete-to-continuous time transformation is +C applied on the system matrices by AB04MD SLICOT subroutine. +C Finally, if requested, the poles and zeros of the system are +C checked. If some of them have positive real parts in the +C continuous-time case (or are not inside the unit disk in the +C complex plane in the discrete-time case), they are exchanged with +C their negatives (or reciprocals, respectively), to preserve the +C frequency response, while getting a minimum phase and stable +C system. This is done by SB10ZP SLICOT subroutine. +C +C REFERENCES +C +C [1] Oppenheim, A.V. and Schafer, R.W. +C Discrete-Time Signal Processing. +C Prentice-Hall Signal Processing Series, 1989. +C +C [2] Balas, G., Doyle, J., Glover, K., Packard, A., and Smith, R. +C Mu-analysis and Synthesis toolbox - User's Guide, +C The Mathworks Inc., Natick, MA, USA, 1998. +C +C CONTRIBUTORS +C +C Asparuh Markovski, Technical University of Sofia, July 2003. +C +C REVISIONS +C +C V. Sima, Research Institute for Informatics, Bucharest, Aug. 2003. +C A. Markovski, Technical University of Sofia, October 2003. +C +C KEYWORDS +C +C Bilinear transformation, frequency response, least-squares +C approximation, stability. +C +C ****************************************************************** +C +C .. Parameters .. + COMPLEX*16 ZZERO, ZONE + PARAMETER ( ZZERO = ( 0.0D+0, 0.0D+0 ), + $ ZONE = ( 1.0D+0, 0.0D+0 ) ) + DOUBLE PRECISION ZERO, ONE, TWO, FOUR, TEN + PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0, + $ FOUR = 4.0D+0, TEN = 1.0D+1 ) + INTEGER HNPTS + PARAMETER ( HNPTS = 2048 ) +C .. +C .. Scalar Arguments .. + INTEGER DISCFL, FLAG, INFO, LDA, LDWORK, LENDAT, + $ LZWORK, N + DOUBLE PRECISION TOL +C .. +C .. Array Arguments .. + INTEGER IWORK(*) + DOUBLE PRECISION A(LDA, *), B(*), C(*), D(*), DWORK(*), + $ IFRDAT(*), OMEGA(*), RFRDAT(*) + COMPLEX*16 ZWORK(*) +C .. +C .. Local Scalars .. + INTEGER CLWMAX, DLWMAX, I, II, INFO2, IP1, IP2, ISTART, + $ ISTOP, IWA0, IWAB, IWBMAT, IWBP, IWBX, IWDME, + $ IWDOMO, IWMAG, IWS, IWVAR, IWXI, IWXR, IWYMAG, + $ K, LW1, LW2, LW3, LW4, MN, N1, N2, P, RANK + DOUBLE PRECISION P1, P2, PI, PW, RAT, TOLB, TOLL + COMPLEX*16 XHAT(HNPTS/2) +C .. +C .. External Functions .. + DOUBLE PRECISION DLAMCH, DLAPY2 + EXTERNAL DLAMCH, DLAPY2 +C .. +C .. External Subroutines .. + EXTERNAL AB04MD, DCOPY, DG01MD, DGELSY, DLASET, DSCAL, + $ SB10ZP, XERBLA +C .. +C .. Intrinsic Functions .. + INTRINSIC ACOS, ATAN, COS, DBLE, DCMPLX, DIMAG, EXP, LOG, + $ MAX, MIN, SIN, SQRT +C +C Test input parameters and workspace. +C + PI = FOUR*ATAN( ONE ) + PW = OMEGA(1) + N1 = N + 1 + N2 = N + N1 +C + INFO = 0 + IF( DISCFL.NE.0 .AND. DISCFL.NE.1 ) THEN + INFO = -1 + ELSE IF( FLAG.NE.0 .AND. FLAG.NE.1 ) THEN + INFO = -2 + ELSE IF ( LENDAT.LT.2 ) THEN + INFO = -3 + ELSE IF ( PW.LT.ZERO ) THEN + INFO = -6 + ELSE IF( N.GT.LENDAT - 1 ) THEN + INFO = -7 + ELSE IF( LDA.LT.MAX( 1, N ) ) THEN + INFO = -9 + ELSE +C + DO 10 K = 2, LENDAT + IF ( OMEGA(K).LT.PW ) + $ INFO = -6 + PW = OMEGA(K) + 10 CONTINUE +C + IF ( DISCFL.EQ.1 .AND. OMEGA(LENDAT).GT.PI ) + $ INFO = -6 + END IF +C + IF ( INFO.EQ.0 ) THEN +C +C Workspace. +C + LW1 = 2*LENDAT + 4*HNPTS + LW2 = LENDAT + 6*HNPTS + MN = MIN( 2*LENDAT, N2 ) +C + IF ( N.GT.0 ) THEN + LW3 = 2*LENDAT*N2 + MAX( 2*LENDAT, N2 ) + + $ MAX( MN + 6*N + 4, 2*MN + 1 ) + ELSE + LW3 = 4*LENDAT + 5 + END IF +C + IF ( FLAG.EQ.0 ) THEN + LW4 = 0 + ELSE + LW4 = MAX( N*N + 5*N, 6*N + 1 + MIN ( 1, N ) ) + END IF +C + DLWMAX = MAX( 2, LW1, LW2, LW3, LW4 ) +C + IF ( N.GT.0 ) THEN + CLWMAX = LENDAT*( N2 + 2 ) + ELSE + CLWMAX = LENDAT + END IF +C + IF ( LDWORK.LT.DLWMAX ) THEN + INFO = -16 + ELSE IF ( LZWORK.LT.CLWMAX ) THEN + INFO = -18 + END IF + END IF +C + IF ( INFO.NE.0 ) THEN +C +C Error return. +C + CALL XERBLA( 'SB10YD', -INFO ) + RETURN + END IF +C +C Set tolerances. +C + TOLB = DLAMCH( 'Epsilon' ) + TOLL = TOL + IF ( TOLL.LE.ZERO ) + $ TOLL = FOUR*DBLE( LENDAT*N )*TOLB +C +C @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ +C +C Workspace usage 1. +C Workspace: need 2*LENDAT + 4*HNPTS. +C + IWDOMO = 1 + IWDME = IWDOMO + LENDAT + IWYMAG = IWDME + 2*HNPTS + IWMAG = IWYMAG + 2*HNPTS +C +C Bilinear transformation. +C + IF ( DISCFL.EQ.0 ) THEN + PW = SQRT( OMEGA(1)*OMEGA(LENDAT) + SQRT( TOLB ) ) +C + DO 20 K = 1, LENDAT + DWORK(IWDME+K-1) = ( OMEGA(K)/PW )**2 + DWORK(IWDOMO+K-1) = + $ ACOS( ( ONE - DWORK(IWDME+K-1) )/ + $ ( ONE + DWORK(IWDME+K-1) ) ) + 20 CONTINUE +C + ELSE + CALL DCOPY( LENDAT, OMEGA, 1, DWORK(IWDOMO), 1 ) + END IF +C +C Linear interpolation. +C + DO 30 K = 1, LENDAT + DWORK(IWMAG+K-1) = DLAPY2( RFRDAT(K), IFRDAT(K) ) + DWORK(IWMAG+K-1) = ( ONE/LOG( TEN ) ) * LOG( DWORK(IWMAG+K-1) ) + 30 CONTINUE +C + DO 40 K = 1, HNPTS + DWORK(IWDME+K-1) = ( K - 1 )*PI/HNPTS + DWORK(IWYMAG+K-1) = ZERO +C + IF ( DWORK(IWDME+K-1).LT.DWORK(IWDOMO) ) THEN + DWORK(IWYMAG+K-1) = DWORK(IWMAG) + ELSE IF ( DWORK(IWDME+K-1).GE.DWORK(IWDOMO+LENDAT-1) ) THEN + DWORK(IWYMAG+K-1) = DWORK(IWMAG+LENDAT-1) + END IF +C + 40 CONTINUE +C + DO 60 I = 2, LENDAT + P1 = HNPTS*DWORK(IWDOMO+I-2)/PI + ONE +C + IP1 = INT( P1 ) + IF ( DBLE( IP1 ).NE.P1 ) + $ IP1 = IP1 + 1 +C + P2 = HNPTS*DWORK(IWDOMO+I-1)/PI + ONE +C + IP2 = INT( P2 ) + IF ( DBLE( IP2 ).NE.P2 ) + $ IP2 = IP2 + 1 +C + DO 50 P = IP1, IP2 - 1 + RAT = DWORK(IWDME+P-1) - DWORK(IWDOMO+I-2) + RAT = RAT/( DWORK(IWDOMO+I-1) - DWORK(IWDOMO+I-2) ) + DWORK(IWYMAG+P-1) = ( ONE - RAT )*DWORK(IWMAG+I-2) + + $ RAT*DWORK(IWMAG+I-1) + 50 CONTINUE +C + 60 CONTINUE +C + DO 70 K = 1, HNPTS + DWORK(IWYMAG+K-1) = EXP( LOG( TEN )*DWORK(IWYMAG+K-1) ) + 70 CONTINUE +C +C Duplicate data around disc. +C + DO 80 K = 1, HNPTS + DWORK(IWDME+HNPTS+K-1) = TWO*PI - DWORK(IWDME+HNPTS-K) + DWORK(IWYMAG+HNPTS+K-1) = DWORK(IWYMAG+HNPTS-K) + 80 CONTINUE +C +C Complex cepstrum to get min phase: +C LOG (Magnitude) +C + DO 90 K = 1, 2*HNPTS + DWORK(IWYMAG+K-1) = TWO*LOG( DWORK(IWYMAG+K-1) ) + 90 CONTINUE +C +C @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ +C +C Workspace usage 2. +C Workspace: need LENDAT + 6*HNPTS. +C + IWXR = IWYMAG + IWXI = IWMAG +C + DO 100 K = 1, 2*HNPTS + DWORK(IWXI+K-1) = ZERO + 100 CONTINUE +C +C IFFT +C + CALL DG01MD( 'I', 2*HNPTS, DWORK(IWXR), DWORK(IWXI), INFO2 ) +C +C Rescale, because DG01MD doesn't do it. +C + CALL DSCAL( HNPTS, ONE/( TWO*HNPTS ), DWORK(IWXR), 1 ) + CALL DSCAL( HNPTS, ONE/( TWO*HNPTS ), DWORK(IWXI), 1 ) +C +C Halve the result at 0. +C + DWORK(IWXR) = DWORK(IWXR)/TWO + DWORK(IWXI) = DWORK(IWXI)/TWO +C +C FFT +C + CALL DG01MD( 'D', HNPTS, DWORK(IWXR), DWORK(IWXI), INFO2 ) +C +C Get the EXP of the result. +C + DO 110 K = 1, HNPTS/2 + XHAT(K) = EXP( DWORK(IWXR+K-1) )* + $ DCMPLX ( COS( DWORK(IWXI+K-1)), SIN( DWORK(IWXI+K-1) ) ) + DWORK(IWDME+K-1) = DWORK(IWDME+2*K-2) + 110 CONTINUE +C +C Interpolate back to original frequency data. +C + ISTART = 1 + ISTOP = LENDAT +C + DO 120 I = 1, LENDAT + ZWORK(I) = ZZERO + IF ( DWORK(IWDOMO+I-1).LE.DWORK(IWDME) ) THEN + ZWORK(I) = XHAT(1) + ISTART = I + 1 + ELSE IF ( DWORK(IWDOMO+I-1).GE.DWORK(IWDME+HNPTS/2-1) ) + $ THEN + ZWORK(I) = XHAT(HNPTS/2) + ISTOP = ISTOP - 1 + END IF + 120 CONTINUE +C + DO 140 I = ISTART, ISTOP + II = HNPTS/2 + 130 CONTINUE + IF ( DWORK(IWDME+II-1).GE.DWORK(IWDOMO+I-1) ) + $ P = II + II = II - 1 + IF ( II.GT.0 ) + $ GOTO 130 + RAT = ( DWORK(IWDOMO+I-1) - DWORK(IWDME+P-2) )/ + $ ( DWORK(IWDME+P-1) - DWORK(IWDME+P-2) ) + ZWORK(I) = RAT*XHAT(P) + ( ONE - RAT )*XHAT(P-1) + 140 CONTINUE +C +C CASE N > 0. +C This is the only allowed case in mu-synthesis subroutines. +C + IF ( N.GT.0 ) THEN +C +C Preparation for frequency identification. +C +C @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ +C +C Complex workspace usage 1. +C Complex workspace: need 2*LENDAT + LENDAT*(N+1). +C + IWA0 = 1 + LENDAT + IWVAR = IWA0 + LENDAT*N1 +C + DO 150 K = 1, LENDAT + IF ( DISCFL.EQ.0 ) THEN + ZWORK(IWVAR+K-1) = DCMPLX( COS( DWORK(IWDOMO+K-1) ), + $ SIN( DWORK(IWDOMO+K-1) ) ) + ELSE + ZWORK(IWVAR+K-1) = DCMPLX( COS( OMEGA(K) ), + $ SIN( OMEGA(K) ) ) + END IF + 150 CONTINUE +C +C Array for DGELSY. +C + DO 160 K = 1, N2 + IWORK(K) = 0 + 160 CONTINUE +C +C Constructing A0. +C + DO 170 K = 1, LENDAT + ZWORK(IWA0+N*LENDAT+K-1) = ZONE + 170 CONTINUE +C + DO 190 I = 1, N + DO 180 K = 1, LENDAT + ZWORK(IWA0+(N-I)*LENDAT+K-1) = + $ ZWORK(IWA0+(N1-I)*LENDAT+K-1)*ZWORK(IWVAR+K-1) + 180 CONTINUE + 190 CONTINUE +C +C @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ +C +C Complex workspace usage 2. +C Complex workspace: need 2*LENDAT + LENDAT*(2*N+1). +C + IWBP = IWVAR + IWAB = IWBP + LENDAT +C +C Constructing BP. +C + DO 200 K = 1, LENDAT + ZWORK(IWBP+K-1) = ZWORK(IWA0+K-1)*ZWORK(K) + 200 CONTINUE +C +C Constructing AB. +C + DO 220 I = 1, N + DO 210 K = 1, LENDAT + ZWORK(IWAB+(I-1)*LENDAT+K-1) = -ZWORK(K)* + $ ZWORK(IWA0+I*LENDAT+K-1) + 210 CONTINUE + 220 CONTINUE +C +C @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ +C +C Workspace usage 3. +C Workspace: need LW3 = 2*LENDAT*(2*N+1) + max(2*LENDAT,2*N+1). +C + IWBX = 1 + 2*LENDAT*N2 + IWS = IWBX + MAX( 2*LENDAT, N2 ) +C +C Constructing AX. +C + DO 240 I = 1, N1 + DO 230 K = 1, LENDAT + DWORK(2*(I-1)*LENDAT+K) = + $ DBLE( ZWORK(IWA0+(I-1)*LENDAT+K-1) ) + DWORK((2*I-1)*LENDAT+K) = + $ DIMAG( ZWORK(IWA0+(I-1)*LENDAT+K-1) ) + 230 CONTINUE + 240 CONTINUE +C + DO 260 I = 1, N + DO 250 K = 1, LENDAT + DWORK(2*N1*LENDAT+2*(I-1)*LENDAT+K) = + $ DBLE( ZWORK(IWAB+(I-1)*LENDAT+K-1) ) + DWORK(2*N1*LENDAT+(2*I-1)*LENDAT+K) = + $ DIMAG( ZWORK(IWAB+(I-1)*LENDAT+K-1) ) + 250 CONTINUE + 260 CONTINUE +C +C Constructing BX. +C + DO 270 K = 1, LENDAT + DWORK(IWBX+K-1) = DBLE( ZWORK(IWBP+K-1) ) + DWORK(IWBX+LENDAT+K-1) = DIMAG( ZWORK(IWBP+K-1) ) + 270 CONTINUE +C +C Estimating X. +C Workspace: need LW3 + max( MN+3*(2*N+1)+1, 2*MN+1 ), +C where MN = min( 2*LENDAT, 2*N+1 ); +C prefer larger. +C + CALL DGELSY( 2*LENDAT, N2, 1, DWORK, 2*LENDAT, DWORK(IWBX), + $ MAX( 2*LENDAT, N2 ), IWORK, TOLL, RANK, + $ DWORK(IWS), LDWORK-IWS+1, INFO2 ) + DLWMAX = MAX( DLWMAX, INT( DWORK(IWS) + IWS - 1 ) ) +C +C Constructing A matrix. +C + DO 280 K = 1, N + A(K,1) = -DWORK(IWBX+N1+K-1) + 280 CONTINUE +C + IF ( N.GT.1 ) + $ CALL DLASET( 'Full', N, N-1, ZERO, ONE, A(1,2), LDA ) +C +C Constructing B matrix. +C + DO 290 K = 1, N + B(K) = DWORK(IWBX+N1+K-1)*DWORK(IWBX) - DWORK(IWBX+K) + 290 CONTINUE +C +C Constructing C matrix. +C + C(1) = -ONE +C + DO 300 K = 2, N + C(K) = ZERO + 300 CONTINUE +C +C Constructing D matrix. +C + D(1) = DWORK(IWBX) +C +C Transform to continuous-time case, if needed. +C Workspace: need max(1,N); +C prefer larger. +C + IF ( DISCFL.EQ.0 ) THEN + CALL AB04MD( 'D', N, 1, 1, ONE, PW, A, LDA, B, LDA, C, 1, + $ D, 1, IWORK, DWORK, LDWORK, INFO2 ) + IF ( INFO2.NE.0 ) THEN + INFO = 1 + RETURN + END IF + DLWMAX = MAX( DLWMAX, INT( DWORK(1) ) ) + END IF +C +C Make all the real parts of the poles and the zeros negative. +C + IF ( FLAG.EQ.1 ) THEN +C +C Workspace: need max(N*N + 5*N, 6*N + 1 + min(1,N)); +C prefer larger. + CALL SB10ZP( DISCFL, N, A, LDA, B, C, D, IWORK, DWORK, + $ LDWORK, INFO ) + IF ( INFO.NE.0 ) + $ RETURN + DLWMAX = MAX( DLWMAX, INT( DWORK(1) ) ) + END IF +C + ELSE +C +C CASE N = 0. +C +C @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ +C +C Workspace usage 4. +C Workspace: need 4*LENDAT. +C + IWBMAT = 1 + 2*LENDAT + IWS = IWBMAT + 2*LENDAT +C +C Constructing AMAT and BMAT. +C + DO 310 K = 1, LENDAT + DWORK(K) = ONE + DWORK(K+LENDAT) = ZERO + DWORK(IWBMAT+K-1) = DBLE( ZWORK(K) ) + DWORK(IWBMAT+LENDAT+K-1) = DIMAG( ZWORK(K) ) + 310 CONTINUE +C +C Estimating D matrix. +C Workspace: need 4*LENDAT + 5; +C prefer larger. +C + IWORK(1) = 0 + CALL DGELSY( 2*LENDAT, 1, 1, DWORK, 2*LENDAT, DWORK(IWBMAT), + $ 2*LENDAT, IWORK, TOLL, RANK, DWORK(IWS), + $ LDWORK-IWS+1, INFO2 ) + DLWMAX = MAX( DLWMAX, INT( DWORK(IWS) + IWS - 1 ) ) +C + D(1) = DWORK(IWBMAT) +C + END IF +C +C @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ +C + DWORK(1) = DLWMAX + DWORK(2) = CLWMAX + RETURN +C +C *** Last line of SB10YD *** + END Property changes on: trunk/octave-forge/extra/control-devel/src/SB10YD.f ___________________________________________________________________ Added: svn:executable + * Added: trunk/octave-forge/extra/control-devel/src/SB10ZP.f =================================================================== --- trunk/octave-forge/extra/control-devel/src/SB10ZP.f (rev 0) +++ trunk/octave-forge/extra/control-devel/src/SB10ZP.f 2011-10-19 16:38:34 UTC (rev 8789) @@ -0,0 +1,339 @@ + SUBROUTINE SB10ZP( DISCFL, N, A, LDA, B, C, D, IWORK, DWORK, + $ LDWORK, INFO ) +C +C SLICOT RELEASE 5.0. +C +C Copyright (c) 2002-2010 NICONET e.V. +C +C This program is free software: you can redistribute it and/or +C modify it under the terms of the GNU General Public License as +C published by the Free Software Foundation, either version 2 of +C the License, or (at your option) any later version. +C +C This program is distributed in the hope that it will be useful, +C but WITHOUT ANY WARRANTY; without even the implied warranty of +C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +C GNU General Public License for more details. +C +C You should have received a copy of the GNU General Public License +C along with this program. If not, see +C <http://www.gnu.org/licenses/>. +C +C PURPOSE +C +C To transform a SISO (single-input single-output) system [A,B;C,D] +C by mirroring its unstable poles and zeros in the boundary of the +C stability domain, thus preserving the frequency response of the +C system, but making it stable and minimum phase. Specifically, for +C a continuous-time system, the positive real parts of its poles +C and zeros are exchanged with their negatives. Discrete-time +C systems are first converted to continuous-time systems using a +C bilinear transformation, and finally converted back. +C +C ARGUMENTS +C +C Input/Output parameters +C +C DISCFL (input) INTEGER +C Indicates the type of the system, as follows: +C = 0: continuous-time system; +C = 1: discrete-time system. +C +C N (input/output) INTEGER +C On entry, the order of the original system. N >= 0. +C On exit, the order of the transformed, minimal system. +C +C A (input/output) DOUBLE PRECISION array, dimension (LDA,N) +C On entry, the leading N-by-N part of this array must +C contain the original system matrix A. +C On exit, the leading N-by-N part of this array contains +C the transformed matrix A, in an upper Hessenberg form. +C +C LDA INTEGER +C The leading dimension of the array A. LDA >= MAX(1,N). +C +C B (input/output) DOUBLE PRECISION array, dimension (N) +C On entry, this array must contain the original system +C vector B. +C On exit, this array contains the transformed vector B. +C +C C (input/output) DOUBLE PRECISION array, dimension (N) +C On entry, this array must contain the original system +C vector C. +C On exit, this array contains the transformed vector C. +C The first N-1 elements are zero (for the exit value of N). +C +C D (input/output) DOUBLE PRECISION array, dimension (1) +C On entry, this array must contain the original system +C scalar D. +C On exit, this array contains the transformed scalar D. +C +C Workspace +C +C IWORK INTEGER array, dimension max(2,N+1) +C +C DWORK DOUBLE PRECISION array, dimension (LDWORK) +C On exit, if INFO = 0, DWORK(1) returns the optimal value +C of LDWORK. +C +C LDWORK INTEGER +C The length of the array DWORK. +C LDWORK >= max(N*N + 5*N, 6*N + 1 + min(1,N)). +C For optimum performance LDWORK should be larger. +C +C Error indicator +C +C INFO INTEGER +C = 0: successful exit; +C < 0: if INFO = -i, the i-th argument had an illegal +C value; +C = 1: if the discrete --> continuous transformation cannot +C be made; +C = 2: if the system poles cannot be found; +C = 3: if the inverse system cannot be found, i.e., D is +C (close to) zero; +C = 4: if the system zeros cannot be found; +C = 5: if the state-space representation of the new +C transfer function T(s) cannot be found; +C = 6: if the continuous --> discrete transformation cannot +C be made. +C +C METHOD +C +C First, if the system is discrete-time, it is transformed to +C continuous-time using alpha = beta = 1 in the bilinear +C transformation implemented in the SLICOT routine AB04MD. +C Then the eigenvalues of A, i.e., the system poles, are found. +C Then, the inverse of the original system is found and its poles, +C i.e., the system zeros, are evaluated. +C The obtained system poles Pi and zeros Zi are checked and if a +C positive real part is detected, it is exchanged by -Pi or -Zi. +C Then the polynomial coefficients of the transfer function +C T(s) = Q(s)/P(s) are found. +C The state-space representation of T(s) is then obtained. +C The system matrices B, C, D are scaled so that the transformed +C system has the same system gain as the original system. +C If the original system is discrete-time, then the result (which is +C continuous-time) is converted back to discrete-time. +C +C CONTRIBUTORS +C +C Asparuh Markovski, Technical University of Sofia, July 2003. +C +C REVISIONS +C +C V. Sima, Research Institute for Informatics, Bucharest, Aug. 2003. +C +C KEYWORDS +C +C Bilinear transformation, stability, state-space representation. +C +C ****************************************************************** +C +C .. Parameters .. + DOUBLE PRECISION ZERO, ONE + PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) +C .. +C .. Scalar Arguments .. + INTEGER DISCFL, INFO, LDA, LDWORK, N +C .. +C .. Array Arguments .. + INTEGER IWORK( * ) + DOUBLE PRECISION A( LDA, * ), B( * ), C( * ), D( * ), DWORK( * ) +C .. +C .. Local Scalars .. + INTEGER I, IDW1, IDW2, IDW3, IMP, IMZ, INFO2, IWA, IWP, + $ IWPS, IWQ, IWQS, LDW1, MAXWRK, REP, REZ + DOUBLE PRECISION RCOND, SCALB, SCALC, SCALD +C .. +C .. Local Arrays .. + INTEGER INDEX(1) +C .. +C .. External Subroutines .. + EXTERNAL AB04MD, AB07ND, DCOPY, DGEEV, DLACPY, DSCAL, + $ MC01PD, TD04AD, XERBLA +C .. +C .. Intrinsic Functions .. + INTRINSIC ABS, INT, MAX, MIN, SIGN, SQRT +C +C Test input parameters and workspace. +C + INFO = 0 + IF ( DISCFL.NE.0 .AND. DISCFL.NE.1 ) THEN + INFO = -1 + ELSE IF ( N.LT.0 ) THEN + INFO = -2 + ELSE IF( LDA.LT.MAX( 1, N ) ) THEN + INFO = -4 + ELSE IF ( LDWORK.LT.MAX( N*N + 5*N, 6*N + 1 + MIN( 1, N ) ) ) THEN + INFO = -10 + END IF +C + IF ( INFO.NE.0 ) THEN +C +C Error return. +C + CALL XERBLA( 'SB10ZP', -INFO ) + RETURN + END IF +C +C Quick return if possible. +C + IF( N.EQ.0 ) THEN + DWORK(1) = ONE + RETURN + END IF +C +C Workspace usage 1. +C + REP = 1 + IMP = REP + N + REZ = IMP + N + IMZ = REZ + N + IWA = REZ + IDW1 = IWA + N*N + LDW1 = LDWORK - IDW1 + 1 +C +C 1. Discrete --> continuous transformation if needed. +C + IF ( DISCFL.EQ.1 ) THEN +C +C Workspace: need max(1,N); +C prefer larger. +C + CALL AB04MD( 'D', N, 1, 1, ONE, ONE, A, LDA, B, LDA, C, 1, + $ D, 1, IWORK, DWORK, LDWORK, INFO2 ) + IF ( INFO2.NE.0 ) THEN + INFO = 1 + RETURN + END IF + MAXWRK = INT( DWORK(1) ) + ELSE + MAXWRK = 0 + END IF +C +C 2. Determine the factors for restoring system gain. +C + SCALD = D(1) + SCALC = SQRT( ABS( SCALD ) ) + SCALB = SIGN( SCALC, SCALD ) +C +C 3. Find the system poles, i.e., the eigenvalues of A. +C Workspace: need N*N + 2*N + 3*N; +C prefer larger. +C + CALL DLACPY( 'Full', N, N, A, LDA, DWORK(IWA), N ) +C + CALL DGEEV( 'N', 'N', N, DWORK(IWA), N, DWORK(REP), DWORK(IMP), + $ DWORK(IDW1), 1, DWORK(IDW1), 1, DWORK(IDW1), LDW1, + $ INFO2 ) + IF ( INFO2.NE.0 ) THEN + INFO = 2 + RETURN + END IF + MAXWRK = MAX( MAXWRK, INT( DWORK(IDW1) + IDW1 - 1 ) ) +C +C 4. Compute the inverse system [Ai, Bi; Ci, Di]. +C Workspace: need N*N + 2*N + 4; +C prefer larger. +C + CALL AB07ND( N, 1, A, LDA, B, LDA, C, 1, D, 1, RCOND, IWORK, + $ DWORK(IDW1), LDW1, INFO2 ) + IF ( INFO2.NE.0 ) THEN + INFO = 3 + RETURN + END IF + MAXWRK = MAX( MAXWRK, INT( DWORK(IDW1) + IDW1 - 1 ) ) +C +C 5. Find the system zeros, i.e., the eigenvalues of Ai. +C Workspace: need 4*N + 3*N; +C prefer larger. +C + IDW1 = IMZ + N + LDW1 = LDWORK - IDW1 + 1 +C + CALL DGEEV( 'N', 'N', N, A, LDA, DWORK(REZ), DWORK(IMZ), + $ DWORK(IDW1), 1, DWORK(IDW1), 1, DWORK(IDW1), LDW1, + $ INFO2 ) + IF ( INFO2.NE.0 ) THEN + INFO = 4 + RETURN + END IF + MAXWRK = MAX( MAXWRK, INT( DWORK(IDW1) + IDW1 - 1 ) ) +C +C 6. Exchange the zeros and the poles with positive real parts with +C their negatives. +C + DO 10 I = 0, N - 1 + IF ( DWORK(REP+I).GT.ZERO ) + $ DWORK(REP+I) = -DWORK(REP+I) + IF ( DWORK(REZ+I).GT.ZERO ) + $ DWORK(REZ+I) = -DWORK(REZ+I) + 10 CONTINUE +C +C Workspace usage 2. +C + IWP = IDW1 + IDW2 = IWP + N + 1 + IWPS = 1 +C +C 7. Construct the nominator and the denominator +C of the system transfer function T( s ) = Q( s )/P( s ). +C 8. Rearrange the coefficients in Q(s) and P(s) because +C MC01PD subroutine produces them in increasing powers of s. +C Workspace: need 6*N + 2. +C + CALL MC01PD( N, DWORK(REP), DWORK(IMP), DWORK(IWP), DWORK(IDW2), + $ INFO2 ) + CALL DCOPY( N+1, DWORK(IWP), -1, DWORK(IWPS), 1 ) +C +C Workspace usage 3. +C + IWQ = IDW1 + IWQS = IWPS + N + 1 + IDW3 = IWQS + N + 1 +C + CALL MC01PD( N, DWORK(REZ), DWORK(IMZ), DWORK(IWQ), DWORK(IDW2), + $ INFO2 ) + CALL DCOPY( N+1, DWORK(IWQ), -1, DWORK(IWQS), 1 ) +C +C 9. Make the conversion T(s) --> [A, B; C, D]. +C Workspace: need 2*N + 2 + N + max(N,3); +C prefer larger. +C + INDEX(1) = N + CALL TD04AD( 'R', 1, 1, INDEX, DWORK(IWPS), 1, DWORK(IWQS), 1, 1, + $ N, A, LDA, B, LDA, C, 1, D, 1, -ONE, IWORK, + $ DWORK(IDW3), LDWORK-IDW3+1, INFO2 ) + IF ( INFO2.NE.0 ) THEN + INFO = 5 + RETURN + END IF + MAXWRK = MAX( MAXWRK, INT( DWORK(IDW3) + IDW3 - 1 ) ) +C +C 10. Scale the transformed system to the previous gain. +C + IF ( N.GT.0 ) THEN + CALL DSCAL( N, SCALB, B, 1 ) + C(N) = SCALC*C(N) + END IF +C + D(1) = SCALD +C +C 11. Continuous --> discrete transformation if needed. +C + IF ( DISCFL.EQ.1 ) THEN + CALL AB04MD( 'C', N, 1, 1, ONE, ONE, A, LDA, B, LDA, C, 1, + $ D, 1, IWORK, DWORK, LDWORK, INFO2 ) + + IF ( INFO2.NE.0 ) THEN + INFO = 6 + RETURN + END IF + END IF +C + DWORK(1) = MAXWRK + RETURN +C +C *** Last line of SB10ZP *** + END Property changes on: trunk/octave-forge/extra/control-devel/src/SB10ZP.f ___________________________________________________________________ Added: svn:executable + * Added: trunk/octave-forge/extra/control-devel/src/TD03AY.f =================================================================== --- trunk/octave-forge/extra/control-devel/src/TD03AY.f (rev 0) +++ trunk/octave-forge/extra/control-devel/src/TD03AY.f 2011-10-19 16:38:34 UTC (rev 8789) @@ -0,0 +1,171 @@ + SUBROUTINE TD03AY( MWORK, PWORK, INDEX, DCOEFF, LDDCOE, UCOEFF, + $ LDUCO1, LDUCO2, N, A, LDA, B, LDB, C, LDC, D, + $ LDD, INFO ) +C +C SLICOT RELEASE 5.0. +C +C Copyright (c) 2002-2010 NICONET e.V. +C +C This program is free software: you can redistribute it and/or +C modify it under the terms of the GNU General Public License as +C published by the Free Software Foundation, either version 2 of +C the License, or (at your option) any later version. +C +C This program is distributed in the hope that it will be useful, +C but WITHOUT ANY WARRANTY; without even the implied warranty of +C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +C GNU General Public License for more details. +C +C You should have received a copy of the GNU General Public License +C along with this program. If not, see +C <http://www.gnu.org/licenses/>. +C +C Calculates a state-space representation for a (PWORK x MWORK) +C transfer matrix given in the form of polynomial row vectors over +C common denominators (not necessarily lcd's). Such a description +C is simply the polynomial matrix representation +C +C T(s) = inv(D(s)) * U(s), +C +C where D(s) is diagonal with (I,I)-th element D:I(s) of degree +C INDEX(I); applying Wolovich's Observable Structure Theorem to +C this left matrix fraction then yields an equivalent state-space +C representation in observable companion form, of order +C N = sum(INDEX(I)). As D(s) is diagonal, the PWORK ordered +C 'non-trivial' columns of C and A are very simply calculated, these +C submatrices being diagonal and (INDEX(I) x 1) - block diagonal, +C respectively: finding B and D is also somewhat simpler than for +C general P(s) as dealt with in TC04AD. Finally, the state-space +C representation obtained here is not necessarily controllable +C (as D(s) and U(s) are not necessarily relatively left prime), but +C it is theoretically completely observable: however, its +C observability matrix may be poorly conditioned, so it is safer +C not to assume observability either. +C +C REVISIONS +C +C May 13, 1998. +C +C KEYWORDS +C +C Coprime matrix fraction, elementary polynomial operations, +C polynomial matrix, state-space representation, transfer matrix. +C +C ****************************************************************** +C + DOUBLE PRECISION ZERO, ONE + PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) +C .. Scalar Arguments .. + INTEGER INFO, LDA, LDB, LDC, LDD, LDDCOE, LDUCO1, + $ LDUCO2, MWORK, N, PWORK +C .. Array Arguments .. + INTEGER INDEX(*) + DOUBLE PRECISION A(LDA,*), B(LDB,*), C(LDC,*), D(LDD,*), + $ DCOEFF(LDDCOE,*), UCOEFF(LDUCO1,LDUCO2,*) +C .. Local Scalars .. + INTEGER I, IA, IBIAS, INDCUR, JA, JMAX1, K + DOUBLE PRECISION ABSDIA, ABSDMX, BIGNUM, DIAG, SMLNUM, UMAX1, + $ TEMP +C .. External Functions .. + INTEGER IDAMAX + DOUBLE PRECISION DLAMCH + EXTERNAL DLAMCH, IDAMAX +C .. External Subroutines .. + EXTERNAL DAXPY, DCOPY, DLASET, DSCAL +C .. Intrinsic Functions .. + INTRINSIC ABS +C .. Executable Statements .. +C + INFO = 0 +C +C Initialize A and C to be zero, apart from 1's on the subdiagonal +C of A. +C + CALL DLASET( 'Upper', N, N, ZERO, ZERO, A, LDA ) + IF ( N.GT.1 ) CALL DLASET( 'Lower', N-1, N-1, ZERO, ONE, A(2,1), + $ LDA ) +C + CALL DLASET( 'Full', PWORK, N, ZERO, ZERO, C, LDC ) +C +C Calculate B and D, as well as 'non-trivial' elements of A and C. +C Check if any leading coefficient of D(s) nearly zero: if so, exit. +C Caution is taken to avoid overflow. +C + SMLNUM = DLAMCH( 'Safe minimum' ) / DLAMCH( 'Precision' ) + BIGNUM = ONE / SMLNUM +C + IBIAS = 2 + JA = 0 +C + DO 20 I = 1, PWORK + ABSDIA = ABS( DCOEFF(I,1) ) + JMAX1 = IDAMAX( MWORK, UCOEFF(I,1,1), LDUCO1 ) + UMAX1 = ABS( UCOEFF(I,JMAX1,1) ) + IF ( ( ABSDIA.LT.SMLNUM ) .OR. + $ ( ABSDIA.LT.ONE .AND. UMAX1.GT.ABSDIA*BIGNUM ) ) THEN +C +C Error return. +C + INFO = I + RETURN + END IF + DIAG = ONE/DCOEFF(I,1) + INDCUR = INDEX(I) + IF ( INDCUR.NE.0 ) THEN + IBIAS = IBIAS + INDCUR + JA = JA + INDCUR + IF ( INDCUR.GE.1 ) THEN + JMAX1 = IDAMAX( INDCUR, DCOEFF(I,2), LDDCOE ) + ABSDMX = ABS( DCOEFF(I,JMAX1) ) + IF ( ABSDIA.GE.ONE ) THEN + IF ( UMAX1.GT.ONE ) THEN + IF ( ( ABSDMX/ABSDIA ).GT.( BIGNUM/UMAX1 ) ) THEN +C +C Error return. +C + INFO = I + RETURN + END IF + END IF + ELSE + IF ( UMAX1.GT.ONE ) THEN + IF ( ABSDMX.GT.( BIGNUM*ABSDIA )/UMAX1 ) THEN +C +C Error return. +C + INFO = I + RETURN + END IF + END IF + END IF + END IF +C +C I-th 'non-trivial' sub-vector of A given from coefficients +C of D:I(s), while I-th row block of B given from this and +C row I of U(s). +C + DO 10 K = 2, INDCUR + 1 + IA = IBIAS - K + TEMP = -DIAG*DCOEFF(I,K) + A(IA,JA) = TEMP +C + CALL DCOPY( MWORK, UCOEFF(I,1,K), LDUCO1, B(IA,1), LDB ) + CALL DAXPY( MWORK, TEMP, UCOEFF(I,1,1), LDUCO1, B(IA,1), + $ LDB ) + 10 CONTINUE +C + IF ( JA.LT.N ) A(JA+1,JA) = ZERO +C +C Finally, I-th 'non-trivial' entry of C and row of D obtained +C also. +C + C(I,JA) = DIAG + END IF +C + CALL DCOPY( MWORK, UCOEFF(I,1,1), LDUCO1, D(I,1), LDD ) + CALL DSCAL( MWORK, DIAG, D(I,1), LDD ) + 20 CONTINUE +C + RETURN +C *** Last line of TD03AY *** + END Property changes on: trunk/octave-forge/extra/control-devel/src/TD03AY.f ___________________________________________________________________ Added: svn:executable + * Added: trunk/octave-forge/extra/control-devel/src/TD04AD.f =================================================================== --- trunk/octave-forge/extra/control-devel/src/TD04AD.f (rev 0) +++ trunk/octave-forge/extra/control-devel/src/TD04AD.f 2011-10-19 16:38:34 UTC (rev 8789) @@ -0,0 +1,425 @@ + SUBROUTINE TD04AD( ROWCOL, M, P, INDEX, DCOEFF, LDDCOE, UCOEFF, + $ LDUCO1, LDUCO2, NR, A, LDA, B, LDB, C, LDC, D, + $ LDD, TOL, IWORK, DWORK, LDWORK, INFO ) +C +C SLICOT RELEASE 5.0. +C +C Copyright (c) 2002-2010 NICONET e.V. +C +C This program is free software: you can redistribute it and/or +C modify it under the terms of the GNU General Public License as +C published by the Free Software Foundation, either version 2 of +C the License, or (at your option) any later version. +C +C This program is distributed in the hope that it will be useful, +C but WITHOUT ANY WARRANTY; without even the implied warranty of +C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +C GNU General Public License for more details. +C +C You should have received a copy of the GNU General Public License +C along with this program. If not, see +C <http://www.gnu.org/licenses/>. +C +C PURPOSE +C +C To find a minimal state-space representation (A,B,C,D) for a +C proper transfer matrix T(s) given as either row or column +C polynomial vectors over denominator polynomials, possibly with +C uncancelled common terms. +C +C ARGUMENTS +C +C Mode Parameters +C +C ROWCOL CHARACTER*1 +C Indicates whether the transfer matrix T(s) is given as +C rows or columns over common denominators as follows: +C = 'R': T(s) is given as rows over common denominators; +C = 'C': T(s) is given as columns over common denominators. +C +C Input/Output Parameters +C +C M (input) INTEGER +C The number of system inputs. M >= 0. +C +C P (input) INTEGER +C The number of system outputs. P >= 0. +C +C INDEX (input) INTEGER array, dimension (porm), where porm = P, +C if ROWCOL = 'R', and porm = M, if ROWCOL = 'C'. +C This array must contain the degrees of the denominator +C polynomials in D(s). +C +C DCOEFF (input) DOUBLE PRECISION array, dimension (LDDCOE,kdcoef), +C where kdcoef = MAX(INDEX(I)) + 1. +C The leading porm-by-kdcoef part of this array must contain +C the coefficients of each denominator polynomial. +C DCOEFF(I,K) is the coefficient in s**(INDEX(I)-K+1) of the +C I-th denominator polynomial in D(s), where +C K = 1,2,...,kdcoef. +C +C LDDCOE INTEGER +C The leading dimension of array DCOEFF. +C LDDCOE >= MAX(1,P) if ROWCOL = 'R'; +C LDDCOE >= MAX(1,M) if ROWCOL = 'C'. +C +C UCOEFF (input) DOUBLE PRECISION array, dimension +C (LDUCO1,LDUCO2,kdcoef) +C The leading P-by-M-by-kdcoef part of this array must +C contain the numerator matrix U(s); if ROWCOL = 'C', this +C array is modified internally but restored on exit, and the +C remainder of the leading MAX(M,P)-by-MAX(M,P)-by-kdcoef +C part is used as internal workspace. +C UCOEFF(I,J,K) is the coefficient in s**(INDEX(iorj)-K+1) +C of polynomial (I,J) of U(s), where K = 1,2,...,kdcoef; +C if ROWCOL = 'R' then iorj = I, otherwise iorj = J. +C Thus for ROWCOL = 'R', U(s) = +C diag(s**INDEX(I))*(UCOEFF(.,.,1)+UCOEFF(.,.,2)/s+...). +C +C LDUCO1 INTEGER +C The leading dimension of array UCOEFF. +C LDUCO1 >= MAX(1,P) if ROWCOL = 'R'; +C LDUCO1 >= MAX(1,M,P) if ROWCOL = 'C'. +C +C LDUCO2 INTEGER +C The second dimension of array UCOEFF. +C LDUCO2 >= MAX(1,M) if ROWCOL = 'R'; +C LDUCO2 >= MAX(1,M,P) if ROWCOL = 'C'. +C +C NR (output) INTEGER +C The order of the resulting minimal realization, i.e. the +C order of th... [truncated message content] |
From: <par...@us...> - 2011-10-19 18:16:33
|
Revision: 8792 http://octave.svn.sourceforge.net/octave/?rev=8792&view=rev Author: paramaniac Date: 2011-10-19 18:16:26 +0000 (Wed, 19 Oct 2011) Log Message: ----------- control-devel: add code for fitfrd (3) Modified Paths: -------------- trunk/octave-forge/extra/control-devel/devel/makefile_ident.m trunk/octave-forge/extra/control-devel/inst/fitfrd.m Added Paths: ----------- trunk/octave-forge/extra/control-devel/devel/test_fitfrd.m trunk/octave-forge/extra/control-devel/src/slsb10yd.cc Modified: trunk/octave-forge/extra/control-devel/devel/makefile_ident.m =================================================================== --- trunk/octave-forge/extra/control-devel/devel/makefile_ident.m 2011-10-19 17:01:28 UTC (rev 8791) +++ trunk/octave-forge/extra/control-devel/devel/makefile_ident.m 2011-10-19 18:16:26 UTC (rev 8792) @@ -24,7 +24,8 @@ MB03UD.f MB01SD.f ## fit state-space model to frequency response data -mkoctfile SB10YD.f DG01MD.f AB04MD.f SB10ZP.f AB07ND.f \ +mkoctfile slsb10yd.cc \ + SB10YD.f DG01MD.f AB04MD.f SB10ZP.f AB07ND.f \ MC01PD.f TD04AD.f TD03AY.f TB01PD.f TB01XD.f \ AB07MD.f TB01UD.f TB01ID.f MB01PD.f MB03OY.f \ MB01QD.f Added: trunk/octave-forge/extra/control-devel/devel/test_fitfrd.m =================================================================== --- trunk/octave-forge/extra/control-devel/devel/test_fitfrd.m (rev 0) +++ trunk/octave-forge/extra/control-devel/devel/test_fitfrd.m 2011-10-19 18:16:26 UTC (rev 8792) @@ -0,0 +1,6 @@ +sys = ss (-1, 1, 1, 0); +n = 1; + +ret0 = fitfrd (sys, n) + +ret1 = fitfrd (sys, n, 1) \ No newline at end of file Modified: trunk/octave-forge/extra/control-devel/inst/fitfrd.m =================================================================== --- trunk/octave-forge/extra/control-devel/inst/fitfrd.m 2011-10-19 17:01:28 UTC (rev 8791) +++ trunk/octave-forge/extra/control-devel/inst/fitfrd.m 2011-10-19 18:16:26 UTC (rev 8792) @@ -16,9 +16,8 @@ ## along with LTI Syncope. If not, see <http://www.gnu.org/licenses/>. ## -*- texinfo -*- -## @deftypefn{Function File} {[@var{K}, @var{N}, @var{gamma}, @var{info}] =} ncfsyn (@var{G}, @var{W1}, @var{W2}, @var{factor}) -## Normalized Coprime Factor (NCF) H-infinity synthesis. -## Compute positive feedback controller using the McFarlane/Glover Loop Shaping Design Procedure. +## @deftypefn{Function File} {[@var{sys}, @var{n}] =} fitfrd (@var{dat}, @var{n}) +## Fit frequency response data with a stable, minimum-phase state-space system. ## ## @strong{Inputs} ## @table @var @@ -68,8 +67,25 @@ ## Created: October 2011 ## Version: 0.1 -function [retsys, n] = fitfrd (sys, n, code) +function [sys, n] = fitfrd (dat, n, flag = 0) + if (nargin == 0 || nargin > 3) + print_usage (); + endif + if (! isa (dat, "frd")) + dat = frd (dat); + endif + + if (! issample (n, 0) || n != round (n)) + error ("fitfrd: second argument must be an integer >= 0"); + endif + [H, w, tsam] = frdata (dat, "vector"); + dt = isdt (dat); + + [a, b, c, d, n] = slsb10yd (real (H), imag (H), w, n, dt, flag); + + sys = ss (a, b, c, d, tsam); + endfunction \ No newline at end of file Added: trunk/octave-forge/extra/control-devel/src/slsb10yd.cc =================================================================== --- trunk/octave-forge/extra/control-devel/src/slsb10yd.cc (rev 0) +++ trunk/octave-forge/extra/control-devel/src/slsb10yd.cc 2011-10-19 18:16:26 UTC (rev 8792) @@ -0,0 +1,159 @@ +/* + +Copyright (C) 2011 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/>. + +Fit FRD with SS model. +Uses SLICOT SB10YD by courtesy of NICONET e.V. +<http://www.slicot.org> + +Author: Lukas Reichlin <luk...@gm...> +Created: October 2011 +Version: 0.1 + +*/ + +#include <octave/oct.h> +#include <f77-fcn.h> +#include <complex> +#include "common.cc" + +extern "C" +{ + int F77_FUNC (sb10yd, SB10YD) + (int& DISCFL, int& FLAG, + int& LENDAT, + double* RFRDAT, double* IFRDAT, + double* OMEGA, + int& N, + double* A, int& LDA, + double* B, + double* C, + double* D, + double& TOL, + int* IWORK, double* DWORK, int& LDWORK, + Complex* ZWORK, int& LZWORK, + int& INFO); +} + +DEFUN_DLD (slsb10yd, args, nargout, + "-*- texinfo -*-\n\ +Slicot SB10YD Release 5.0\n\ +No argument checking.\n\ +For internal use only.") +{ + int nargin = args.length (); + octave_value_list retval; + + if (nargin != 6) + { + print_usage (); + } + else + { + // arguments in + Matrix rfrdat = args(0).matrix_value (); + Matrix ifrdat = args(1).matrix_value (); + Matrix omega = args(2).matrix_value (); + int n = args(3).int_value (); + int discfl = args(4).int_value (); + int flag = args(5).int_value (); + + int lendat = omega.rows (); // number of frequencies + + int lda = max (1, n); + + // arguments out + Matrix a (lda, n); + Matrix b (n, 1); + Matrix c (1, n); + Matrix d (1, 1); + + // workspace + int liwork = max (2, 2*n + 1); + int ldwork; + int lzwork; + int hnpts = 2048; + + int lw1 = 2*lendat + 4*hnpts; + int lw2 = lendat + 6*hnpts; + int mn = min (2*lendat, 2*n+1); + int lw3; + int lw4; + + if (n > 0) + { + lw3 = 2*lendat*(2*n+1) + max (2*lendat, 2*n+1) + + max (mn + 6*n + 4, 2*mn + 1); + lzwork = lendat*(2*n+3); + } + else + { + lw3 = 4*lendat + 5; + lzwork = lendat; + } + + if (flag == 1) + lw4 = max (n*n + 5*n, 6*n + 1 + min (1, n)); + else + lw4 = 0; + + ldwork = max (2, lw1, lw2, lw3, lw4); + + OCTAVE_LOCAL_BUFFER (int, iwork, liwork); + OCTAVE_LOCAL_BUFFER (double, dwork, ldwork); + OCTAVE_LOCAL_BUFFER (Complex, zwork, lzwork); + + // tolerance + double tol = 0; + + // error indicator + int info; + + + // SLICOT routine SB10YD + F77_XFCN (sb10yd, SB10YD, + (discfl, flag, + lendat, + rfrdat.fortran_vec (), ifrdat.fortran_vec (), + omega.fortran_vec (), + n, + a.fortran_vec (), lda, + b.fortran_vec (), + c.fortran_vec (), + d.fortran_vec (), + tol, + iwork, dwork, ldwork, + zwork, lzwork, + info)); + + if (f77_exception_encountered) + error ("fitfrd: slsb10yd: exception in SLICOT subroutine SB10YD"); + + if (info != 0) + error ("fitfrd: slsb10yd: SB10YD returned info = %d", info); + + // return values + retval(0) = a; + retval(1) = b; + retval(2) = c; + retval(3) = d; + retval(4) = octave_value (n); + } + + return retval; +} This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2011-10-19 18:33:09
|
Revision: 8793 http://octave.svn.sourceforge.net/octave/?rev=8793&view=rev Author: paramaniac Date: 2011-10-19 18:33:03 +0000 (Wed, 19 Oct 2011) Log Message: ----------- control-devel: touch up fitfrd Modified Paths: -------------- trunk/octave-forge/extra/control-devel/inst/fitfrd.m trunk/octave-forge/extra/control-devel/src/Makefile Modified: trunk/octave-forge/extra/control-devel/inst/fitfrd.m =================================================================== --- trunk/octave-forge/extra/control-devel/inst/fitfrd.m 2011-10-19 18:16:26 UTC (rev 8792) +++ trunk/octave-forge/extra/control-devel/inst/fitfrd.m 2011-10-19 18:33:03 UTC (rev 8793) @@ -17,45 +17,29 @@ ## -*- texinfo -*- ## @deftypefn{Function File} {[@var{sys}, @var{n}] =} fitfrd (@var{dat}, @var{n}) +## @deftypefnx{Function File} {[@var{sys}, @var{n}] =} fitfrd (@var{dat}, @var{n}, @var{flag}) ## Fit frequency response data with a stable, minimum-phase state-space system. ## ## @strong{Inputs} ## @table @var -## @item G -## LTI model of plant. -## @item W1 -## LTI model of precompensator. Model must be SISO or of appropriate size. -## An identity matrix is taken if @var{W1} is not specified or if an empty model -## @code{[]} is passed. -## @item W2 -## LTI model of postcompensator. Model must be SISO or of appropriate size. -## An identity matrix is taken if @var{W2} is not specified or if an empty model -## @code{[]} is passed. -## @item factor -## @code{factor = 1} implies that an optimal controller is required. -## @code{factor > 1} implies that a suboptimal controller is required, -## achieving a performance that is @var{factor} times less than optimal. -## Default value is 1. +## @item dat +## LTI model containing frequency response data of a SISO system. +## @item n +## The desired order of the system to be fitted. @code{n <= length(dat.w)}. +## @item flag = 0 +## The system zeros and poles are not constrained. Default value +## @item flag = 1 +## The system zeros and poles will have negative real parts in the +## continuous-time case, or moduli less than 1 in the discrete-time case. ## @end table ## ## @strong{Outputs} ## @table @var -## @item K +## @item sys ## State-space model of the H-infinity loop-shaping controller. -## @item N -## State-space model of the closed loop depicted below. -## @item gamma -## L-infinity norm of @var{N}. -## @item info -## Structure containing additional information. -## @item info.emax -## Nugap robustness. @code{emax = inv (gamma)}. -## @item info.Gs -## Shaped plant. @code{Gs = W2 * G * W1}. -## @item info.Ks -## Controller for shaped plant. @code{Ks = ncfsyn (Gs)}. -## @item info.rcond -## Estimates of the reciprocal condition numbers of the Riccati equations. +## @item n +## The order of the obtained system. The value of @var{n} +## could only be modified if inputs @code{n > 0} and @code{flag = 1}. ## @end table ## ## @strong{Algorithm}@* @@ -84,8 +68,12 @@ [H, w, tsam] = frdata (dat, "vector"); dt = isdt (dat); - [a, b, c, d, n] = slsb10yd (real (H), imag (H), w, n, dt, flag); + if (n > length (w)) + error ("fitfrd: require n <= length (dat.w)"); + endif + [a, b, c, d, n] = slsb10yd (real (H), imag (H), w, n, dt, logical (flag)); + sys = ss (a, b, c, d, tsam); endfunction \ No newline at end of file Modified: trunk/octave-forge/extra/control-devel/src/Makefile =================================================================== --- trunk/octave-forge/extra/control-devel/src/Makefile 2011-10-19 18:16:26 UTC (rev 8792) +++ trunk/octave-forge/extra/control-devel/src/Makefile 2011-10-19 18:33:03 UTC (rev 8793) @@ -1,4 +1,4 @@ -all: slab09jd.oct +all: slab09jd.oct slsb10yd.oct # TODO: leading and trailing underscores for sl* functions # (__sl*__.oct) would be nice, but this can be an issue @@ -16,5 +16,13 @@ MB04ND.f MB04OD.f SB03OR.f SB03OY.f MB04NY.f \ MB04OY.f SB03OV.f +# fit state-space model to frequency response data +slsb10yd.oct: slsb10yd.cc + mkoctfile slsb10yd.cc \ + SB10YD.f DG01MD.f AB04MD.f SB10ZP.f AB07ND.f \ + MC01PD.f TD04AD.f TD03AY.f TB01PD.f TB01XD.f \ + AB07MD.f TB01UD.f TB01ID.f MB01PD.f MB03OY.f \ + MB01QD.f + clean: rm *.o core octave-core *.oct *~ This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2011-10-20 19:50:12
|
Revision: 8800 http://octave.svn.sourceforge.net/octave/?rev=8800&view=rev Author: paramaniac Date: 2011-10-20 19:50:05 +0000 (Thu, 20 Oct 2011) Log Message: ----------- control-devel: add/update draft code Modified Paths: -------------- trunk/octave-forge/extra/control-devel/devel/test_ab09jd.m trunk/octave-forge/extra/control-devel/src/slab09jd.cc Added Paths: ----------- trunk/octave-forge/extra/control-devel/inst/hnamodred.m Modified: trunk/octave-forge/extra/control-devel/devel/test_ab09jd.m =================================================================== --- trunk/octave-forge/extra/control-devel/devel/test_ab09jd.m 2011-10-20 16:56:06 UTC (rev 8799) +++ trunk/octave-forge/extra/control-devel/devel/test_ab09jd.m 2011-10-20 19:50:05 UTC (rev 8800) @@ -29,10 +29,34 @@ dv = [ 1 ]; +[ar, br, cr, dr] = slab09jd (a, b, c, d, 0, 0, 0, 1, 0.0, \ + 1, av, bv, cv, dv, \ + 0, [], [], [], [], \ + 2, 1e-1, 1e-14) + + +%{ + 0, 0.0, \ + 1, 0, 2, 0, 0, 1, \ + 1e-1, 1e-14) + [ar, br, cr, dr, nr] = slab09jd (a, b, c, d, dt, scaled, nr, ordsel, alpha, \ + jobv, av, bv, cv, dv, \ + jobw, aw, bw, cw, dw, \ + jobinv, tol1, tol2); +%} +%{ +sys = ss (a, b, c, d); +sysv = ss (av, bv, cv, dv); + +sysr = hnamodred (sys, 0, sysv, []) +%} + +%{ [ar, br, cr, dr] = slab09jd (a, b, c, d, av, bv, cv, dv, [], [], [], [], 0, 0.0, \ 1, 0, 2, 0, 0, 1, \ 1e-1, 1e-14) +%} %{ The reduced state dynamics matrix Ar is Added: trunk/octave-forge/extra/control-devel/inst/hnamodred.m =================================================================== --- trunk/octave-forge/extra/control-devel/inst/hnamodred.m (rev 0) +++ trunk/octave-forge/extra/control-devel/inst/hnamodred.m 2011-10-20 19:50:05 UTC (rev 8800) @@ -0,0 +1,94 @@ +function [sysr, nr] = hnamodred (sys, nr = 0, sysv = [], sysw = []) + + [a, b, c, d, tsam, scaled] = ssdata (sys); + dt = isdt (sys); + + if (isempty (sysv)) + av = bv = cv = dv = []; + jobv = 0; + else + sysv = ss (sysv); + [av, bv, cv, dv] = ssdata (sysv); + jobv = 1; + endif + + if (isempty (sysw)) + aw = bw = cw = dw = []; + jobw = 0; + else + sysw = ss (sysw); + [aw, bw, cw, dw] = ssdata (sysw); + jobw = 1; + endif + + jobinv = 2; + tol1 = 1e-1; + tol2 = 1e-14; + alpha = 0.0; + ordsel = 1; + + [ar, br, cr, dr, nr] = slab09jd (a, b, c, d, dt, scaled, nr, ordsel, alpha, \ + jobv, av, bv, cv, dv, \ + jobw, aw, bw, cw, dw, \ + jobinv, tol1, tol2); + + sysr = ss (ar, br, cr, dr, tsam); + +endfunction + + +%!shared Mo, Me +%! a = [ -3.8637 -7.4641 -9.1416 -7.4641 -3.8637 -1.0000 +%! 1.0000, 0 0 0 0 0 +%! 0 1.0000 0 0 0 0 +%! 0 0 1.0000 0 0 0 +%! 0 0 0 1.0000 0 0 +%! 0 0 0 0 1.0000 0 ]; +%! +%! b = [ 1 +%! 0 +%! 0 +%! 0 +%! 0 +%! 0 ]; +%! +%! c = [ 0 0 0 0 0 1 ]; +%! +%! d = [ 0 ]; +%! +%! sys = ss (a, b, c, d); +%! +%! av = [ 0.2000 -1.0000 +%! 1.0000 0 ]; +%! +%! bv = [ 1 +%! 0 ]; +%! +%! cv = [ -1.8000 0 ]; +%! +%! dv = [ 1 ]; +%! +%! sysv = ss (av, bv, cv, dv); +%! +%! sysr = hnamodred (sys, 0, sysv, []); +%! [ao, bo, co, do] = ssdata (sysr); +%! +%! ae = [ -0.2391 0.3072 1.1630 1.1967 +%! -2.9709 -0.2391 2.6270 3.1027 +%! 0.0000 0.0000 -0.5137 -1.2842 +%! 0.0000 0.0000 0.1519 -0.5137 ]; +%! +%! be = [ -1.0497 +%! -3.7052 +%! 0.8223 +%! 0.7435 ]; +%! +%! ce = [ -0.4466 0.0143 -0.4780 -0.2013 ]; +%! +%! de = [ 0.0219 ]; +%! +%! Mo = [ao, bo; co, do]; +%! Me = [ae, be; ce, de]; +%! +%!assert (Mo, Me, 1e-4); + Modified: trunk/octave-forge/extra/control-devel/src/slab09jd.cc =================================================================== --- trunk/octave-forge/extra/control-devel/src/slab09jd.cc 2011-10-20 16:56:06 UTC (rev 8799) +++ trunk/octave-forge/extra/control-devel/src/slab09jd.cc 2011-10-20 19:50:05 UTC (rev 8800) @@ -87,26 +87,25 @@ Matrix c = args(2).matrix_value (); Matrix d = args(3).matrix_value (); - Matrix av = args(4).matrix_value (); - Matrix bv = args(5).matrix_value (); - Matrix cv = args(6).matrix_value (); - Matrix dv = args(7).matrix_value (); - - Matrix aw = args(8).matrix_value (); - Matrix bw = args(9).matrix_value (); - Matrix cw = args(10).matrix_value (); - Matrix dw = args(11).matrix_value (); + const int idico = args(4).int_value (); + const int iequil = args(5).int_value (); + int nr = args(6).int_value (); + const int iordsel = args(7).int_value (); + double alpha = args(8).double_value (); + + const int ijobv = args(9).int_value (); + Matrix av = args(10).matrix_value (); + Matrix bv = args(11).matrix_value (); + Matrix cv = args(12).matrix_value (); + Matrix dv = args(13).matrix_value (); - int nr = args(12).int_value (); - double alpha = args(13).double_value (); - - const int ijobv = args(14).int_value (); - const int ijobw = args(15).int_value (); - const int ijobinv = args(16).int_value (); - const int idico = args(17).int_value (); - const int iequil = args(18).int_value (); - const int iordsel = args(19).int_value (); - + const int ijobw = args(14).int_value (); + Matrix aw = args(15).matrix_value (); + Matrix bw = args(16).matrix_value (); + Matrix cw = args(17).matrix_value (); + Matrix dw = args(18).matrix_value (); + + const int ijobinv = args(19).int_value (); double tol1 = args(20).double_value (); double tol2 = args(21).double_value (); @@ -308,6 +307,7 @@ retval(1) = b; retval(2) = c; retval(3) = d; + retval(4) = octave_value (nr); // retval(0) = hsv; // retval(1) = octave_value (ns); } This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2011-10-25 21:30:26
|
Revision: 8862 http://octave.svn.sourceforge.net/octave/?rev=8862&view=rev Author: paramaniac Date: 2011-10-25 21:30:19 +0000 (Tue, 25 Oct 2011) Log Message: ----------- control-devel: compare different results Modified Paths: -------------- trunk/octave-forge/extra/control-devel/inst/hnamodred.m Added Paths: ----------- trunk/octave-forge/extra/control-devel/devel/compare_results_hnamodred.m Added: trunk/octave-forge/extra/control-devel/devel/compare_results_hnamodred.m =================================================================== --- trunk/octave-forge/extra/control-devel/devel/compare_results_hnamodred.m (rev 0) +++ trunk/octave-forge/extra/control-devel/devel/compare_results_hnamodred.m 2011-10-25 21:30:19 UTC (rev 8862) @@ -0,0 +1,33 @@ + Mo = [ + + -0.23915 0.30723 1.16297 -1.19671 1.04965 + -2.97091 -0.23915 2.62702 -3.10273 3.70515 + 0.00000 0.00000 -0.51368 1.28421 -0.82227 + 0.00000 0.00000 -0.15189 -0.51368 0.74348 + 0.44660 -0.01427 0.47803 -0.20129 0.02190 +]; + + Me = [ + + -0.23910 0.30720 1.16300 1.19670 -1.04970 + -2.97090 -0.23910 2.62700 3.10270 -3.70520 + 0.00000 0.00000 -0.51370 -1.28420 0.82230 + 0.00000 0.00000 0.15190 -0.51370 0.74350 + -0.44660 0.01430 -0.47800 -0.20130 0.02190 +]; + +syso = ss (Mo(1:4, 1:4), Mo(1:4, 5), Mo(5, 1:4), Mo(5, 5)) + +syse = ss (Me(1:4, 1:4), Me(1:4, 5), Me(5, 1:4), Me(5, 5)) + +figure (1) +bode (syso) + +figure (2) +bode (syse) + +figure (3) +step (syso) + +figure (4) +step (syse) \ No newline at end of file Modified: trunk/octave-forge/extra/control-devel/inst/hnamodred.m =================================================================== --- trunk/octave-forge/extra/control-devel/inst/hnamodred.m 2011-10-25 15:20:09 UTC (rev 8861) +++ trunk/octave-forge/extra/control-devel/inst/hnamodred.m 2011-10-25 21:30:19 UTC (rev 8862) @@ -133,6 +133,8 @@ endswitch endfor + ## TODO: handle jobv, jobw, (jobinv) + ## perform model order reduction [ar, br, cr, dr, nr] = slab09jd (a, b, c, d, dt, scaled, nr, ordsel, alpha, \ jobv, av, bv, cv, dv, \ This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2011-10-26 18:07:16
|
Revision: 8865 http://octave.svn.sourceforge.net/octave/?rev=8865&view=rev Author: paramaniac Date: 2011-10-26 18:07:10 +0000 (Wed, 26 Oct 2011) Log Message: ----------- control-devel: touch up hnamodred Modified Paths: -------------- trunk/octave-forge/extra/control-devel/inst/hnamodred.m trunk/octave-forge/extra/control-devel/src/slab09jd.cc Modified: trunk/octave-forge/extra/control-devel/inst/hnamodred.m =================================================================== --- trunk/octave-forge/extra/control-devel/inst/hnamodred.m 2011-10-26 16:23:52 UTC (rev 8864) +++ trunk/octave-forge/extra/control-devel/inst/hnamodred.m 2011-10-26 18:07:10 UTC (rev 8865) @@ -178,7 +178,7 @@ %! %! DV = [ 1 ]; %! -%! sysv = ss (AV, BV, CV, DV); +%! sysv = ss (AV, BV, CV, DV, "scaled", true); %! %! sysr = hnamodred (sys, "left", sysv, "tol1", 1e-1, "tol2", 1e-14); %! [Ao, Bo, Co, Do] = ssdata (sysr); Modified: trunk/octave-forge/extra/control-devel/src/slab09jd.cc =================================================================== --- trunk/octave-forge/extra/control-devel/src/slab09jd.cc 2011-10-26 16:23:52 UTC (rev 8864) +++ trunk/octave-forge/extra/control-devel/src/slab09jd.cc 2011-10-26 18:07:10 UTC (rev 8865) @@ -380,6 +380,29 @@ } } } + + if (iwarn != 0) + { + switch (iwarn) + { + case 1: + warning ("hnamodred: 1: with ORDSEL = 'F', the selected order NR is greater\ + than NSMIN, the sum of the order of the\ + ALPHA-unstable part and the order of a minimal\ + realization of the ALPHA-stable part of the given\ + system. In this case, the resulting NR is set equal\ + to NSMIN."); + break; + case 2: + warning ("hnamodred: 2: with ORDSEL = 'F', the selected order NR is less\ + than the order of the ALPHA-unstable part of the\ + given system. In this case NR is set equal to the\ + order of the ALPHA-unstable part."); + break; + default: + warning ("hnamodred: unknown warning, iwarn = %d", info); + } + } // resize a.resize (nr, nr); This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2011-10-26 18:32:59
|
Revision: 8866 http://octave.svn.sourceforge.net/octave/?rev=8866&view=rev Author: paramaniac Date: 2011-10-26 18:32:52 +0000 (Wed, 26 Oct 2011) Log Message: ----------- control-devel: add oct-file for bstmodred Modified Paths: -------------- trunk/octave-forge/extra/control-devel/devel/makefile_modred.m Added Paths: ----------- trunk/octave-forge/extra/control-devel/src/slab09hd.cc Modified: trunk/octave-forge/extra/control-devel/devel/makefile_modred.m =================================================================== --- trunk/octave-forge/extra/control-devel/devel/makefile_modred.m 2011-10-26 18:07:10 UTC (rev 8865) +++ trunk/octave-forge/extra/control-devel/devel/makefile_modred.m 2011-10-26 18:32:52 UTC (rev 8866) @@ -23,7 +23,9 @@ MB04ND.f MB04OD.f SB03OR.f SB03OY.f MB04NY.f \ MB04OY.f SB03OV.f -mkoctfile AB09HD.f TB01ID.f AB04MD.f TB01KD.f AB09HY.f \ +mkoctfile "-Wl,-framework" "-Wl,vecLib" \ + slab09hd.cc \ + AB09HD.f TB01ID.f AB04MD.f TB01KD.f AB09HY.f \ AB09IX.f MB03UD.f SB02MD.f AB09DD.f TB01LD.f \ SB03OU.f MA02AD.f MB03QX.f select.f SB03OT.f \ SB02MR.f SB02MS.f MB03QD.f SB02MU.f SB02MV.f \ Added: trunk/octave-forge/extra/control-devel/src/slab09hd.cc =================================================================== --- trunk/octave-forge/extra/control-devel/src/slab09hd.cc (rev 0) +++ trunk/octave-forge/extra/control-devel/src/slab09hd.cc 2011-10-26 18:32:52 UTC (rev 8866) @@ -0,0 +1,271 @@ +/* + +Copyright (C) 2011 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/>. + +Model reduction based on balanced stochastic truncation method. +Uses SLICOT AB09HD by courtesy of NICONET e.V. +<http://www.slicot.org> + +Author: Lukas Reichlin <luk...@gm...> +Created: October 2011 +Version: 0.1 + +*/ + +#include <octave/oct.h> +#include <f77-fcn.h> +#include "common.cc" + +extern "C" +{ + int F77_FUNC (ab09hd, AB09HD) + (char& DICO, char& JOB, char& EQUIL, char& ORDSEL, + int& N, int& M, int& P, + int& NR, + double& ALPHA, double& BETA, + double* A, int& LDA, + double* B, int& LDB, + double* C, int& LDC, + double* D, int& LDD, + int& NS, + double* HSV, + double& TOL1, double& TOL2, + int* IWORK, + double* DWORK, int& LDWORK, + bool* BWORK, + int& IWARN, int& INFO); +} + +DEFUN_DLD (slab09hd, args, nargout, + "-*- texinfo -*-\n\ +Slicot AB09HD Release 5.0\n\ +No argument checking.\n\ +For internal use only.") +{ + int nargin = args.length (); + octave_value_list retval; + + if (nargin != 13) + { + print_usage (); + } + else + { + // arguments in + char dico; + char job; + char equil; + char ordsel; + + Matrix a = args(0).matrix_value (); + Matrix b = args(1).matrix_value (); + Matrix c = args(2).matrix_value (); + Matrix d = args(3).matrix_value (); + + const int idico = args(4).int_value (); + const int iequil = args(5).int_value (); + const int ijob = args(6).int_value (); + + int nr = args(7).int_value (); + const int iordsel = args(8).int_value (); + + double alpha = args(9).double_value (); + double beta = args(10).double_value (); + double tol1 = args(11).double_value (); + double tol2 = args(12).double_value (); + + switch (ijob) + { + case 0: + job = 'B'; + break; + case 1: + job = 'F'; + break; + case 2: + job = 'S'; + break; + case 3: + job = 'P'; + break; + default: + error ("slab09hd: argument job invalid"); + } + + if (idico == 0) + dico = 'C'; + else + dico = 'D'; + + if (iequil == 0) + equil = 'S'; + else + equil = 'N'; + + if (iordsel == 0) + ordsel = 'F'; + else + ordsel = 'A'; + + int n = a.rows (); // n: number of states + int m = b.columns (); // m: number of inputs + int p = c.rows (); // p: number of outputs + + int lda = max (1, n); + int ldb = max (1, n); + int ldc = max (1, p); + int ldd = max (1, p); + + // arguments out + int ns; + ColumnVector hsv (n); + + // workspace + int liwork = max (1, 2*n); + int mb; + + if (beta == 0) + mb = m; + else + mb = m + p; + + int ldwork = 2*n*n + mb*(n+p) + + max (2, n*(max (n,mb,p)+5), 2*n*p + max (p*(mb+2), 10*n*(n+1))); + + OCTAVE_LOCAL_BUFFER (int, iwork, liwork); + OCTAVE_LOCAL_BUFFER (double, dwork, ldwork); + OCTAVE_LOCAL_BUFFER (bool, bwork, 2*n); + + // error indicators + int iwarn = 0; + int info = 0; + + + // SLICOT routine AB09HD + F77_XFCN (ab09hd, AB09HD, + (dico, job, equil, ordsel, + n, m, p, + nr, + alpha, beta, + a.fortran_vec (), lda, + b.fortran_vec (), ldb, + c.fortran_vec (), ldc, + d.fortran_vec (), ldd, + ns, + hsv.fortran_vec (), + tol1, tol2, + iwork, + dwork, ldwork, + bwork, + iwarn, info)); + + if (f77_exception_encountered) + error ("bstmodred: exception in SLICOT subroutine AB09HD"); + + if (info != 0) + { + if (info < 0) + error ("bstmodred: the %d-th argument had an invalid value", info); + else + { + switch (info) + { + // FIXME: The code below looks nice, but the error message does not + // because there is much white space after each line break + case 1: + error ("bstmodred: 1: the computation of the ordered real Schur form of A\ + failed"); + case 2: + error ("bstmodred: 2: the reduction of the Hamiltonian matrix to real\ + Schur form failed"); + case 3: + error ("bstmodred: 3: the reordering of the real Schur form of the\ + Hamiltonian matrix failed"); + case 4: + error ("bstmodred: 4: the Hamiltonian matrix has less than N stable\ + eigenvalues"); + case 5: + error ("bstmodred: 5: the coefficient matrix U11 in the linear system\ + X*U11 = U21 to determine X is singular to working\ + precision"); + case 6: + error ("bstmodred: 6: BETA = 0 and D has not a maximal row rank"); + case 7: + error ("bstmodred: 7: the computation of Hankel singular values failed"); + case 8: + error ("bstmodred: 8: the separation of the ALPHA-stable/unstable diagonal\ + blocks failed because of very close eigenvalues"); + case 9: + error ("bstmodred: 9: the resulting order of reduced stable part is less\ + than the number of unstable zeros of the stable\ + part"); + default: + error ("bstmodred: unknown error, info = %d", info); + } + } + } + + if (iwarn != 0) + { + switch (iwarn) + { + case 1: + warning ("bstmodred: 1: with ORDSEL = 'F', the selected order NR is greater\ + than NSMIN, the sum of the order of the\ + ALPHA-unstable part and the order of a minimal\ + realization of the ALPHA-stable part of the given\ + system; in this case, the resulting NR is set equal\ + to NSMIN."); + break; + case 2: + warning ("bstmodred: 2: with ORDSEL = 'F', the selected order NR corresponds\ + to repeated singular values for the ALPHA-stable\ + part, which are neither all included nor all\ + excluded from the reduced model; in this case, the\ + resulting NR is automatically decreased to exclude\ + all repeated singular values."); + break; + case 3: + warning ("bstmodred: 3: with ORDSEL = 'F', the selected order NR is less\ + than the order of the ALPHA-unstable part of the\ + given system; in this case NR is set equal to the\ + order of the ALPHA-unstable part."); + break; + default: + warning ("bstmodred: unknown warning, iwarn = %d", info); + } + } + + // resize + a.resize (nr, nr); + b.resize (nr, m); + c.resize (p, nr); + hsv.resize (ns); + + // return values + retval(0) = a; + retval(1) = b; + retval(2) = c; + retval(3) = d; + retval(4) = octave_value (nr); + // retval(0) = hsv; + // retval(1) = octave_value (ns); + } + + return retval; +} This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2011-10-26 19:13:12
|
Revision: 8868 http://octave.svn.sourceforge.net/octave/?rev=8868&view=rev Author: paramaniac Date: 2011-10-26 19:13:06 +0000 (Wed, 26 Oct 2011) Log Message: ----------- control-devel: add an early version of bstmodred Modified Paths: -------------- trunk/octave-forge/extra/control-devel/inst/hnamodred.m Added Paths: ----------- trunk/octave-forge/extra/control-devel/devel/test_bstmodred.m trunk/octave-forge/extra/control-devel/inst/bstmodred.m Added: trunk/octave-forge/extra/control-devel/devel/test_bstmodred.m =================================================================== --- trunk/octave-forge/extra/control-devel/devel/test_bstmodred.m (rev 0) +++ trunk/octave-forge/extra/control-devel/devel/test_bstmodred.m 2011-10-26 19:13:06 UTC (rev 8868) @@ -0,0 +1,28 @@ +A = [ -0.04165 0.0000 4.9200 -4.9200 0.0000 0.0000 0.0000 + -5.2100 -12.500 0.0000 0.0000 0.0000 0.0000 0.0000 + 0.0000 3.3300 -3.3300 0.0000 0.0000 0.0000 0.0000 + 0.5450 0.0000 0.0000 0.0000 -0.5450 0.0000 0.0000 + 0.0000 0.0000 0.0000 4.9200 -0.04165 0.0000 4.9200 + 0.0000 0.0000 0.0000 0.0000 -5.2100 -12.500 0.0000 + 0.0000 0.0000 0.0000 0.0000 0.0000 3.3300 -3.3300 ]; + +B = [ 0.0000 0.0000 + 12.500 0.0000 + 0.0000 0.0000 + 0.0000 0.0000 + 0.0000 0.0000 + 0.0000 12.500 + 0.0000 0.0000 ]; + +C = [ 1.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 + 0.0000 0.0000 0.0000 1.0000 0.0000 0.0000 0.0000 + 0.0000 0.0000 0.0000 0.0000 1.0000 0.0000 0.0000 ]; + +D = [ 0.0000 0.0000 + 0.0000 0.0000 + 0.0000 0.0000 ]; + +sys = ss (A, B, C, D, "scaled", true); + +sysr = bstmodred (sys, "beta", 1.0, "tol1", 0.1, "tol2", 0.0) +[Ao, Bo, Co, Do] = ssdata (sysr); Added: trunk/octave-forge/extra/control-devel/inst/bstmodred.m =================================================================== --- trunk/octave-forge/extra/control-devel/inst/bstmodred.m (rev 0) +++ trunk/octave-forge/extra/control-devel/inst/bstmodred.m 2011-10-26 19:13:06 UTC (rev 8868) @@ -0,0 +1,195 @@ +## Copyright (C) 2011 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{sysr}, @var{nr}] =} hnamodred (@var{sys}, @dots{}) +## Model order reduction by frequency weighted optimal Hankel-norm approximation method. +## +## @strong{Inputs} +## @table @var +## @item sys +## LTI model to be reduced. +## @item @dots{} +## Pairs of properties and values. +## TODO: describe options. +## @end table +## +## @strong{Outputs} +## @table @var +## @item sysr +## Reduced order state-space model. +## @item nr +## The order of the obtained system @var{sysr}. +## @end table +## +## @strong{Algorithm}@* +## Uses SLICOT AB09HD by courtesy of +## @uref{http://www.slicot.org, NICONET e.V.} +## @end deftypefn + +## Author: Lukas Reichlin <luk...@gm...> +## Created: October 2011 +## Version: 0.1 + +function [sysr, nr] = bstmodred (sys, varargin) + + if (nargin == 0) + print_usage (); + endif + + if (! isa (sys, "lti")) + error ("bstmodred: first argument must be an LTI system"); + endif + + if (rem (nargin-1, 2)) + error ("bstmodred: properties and values must come in pairs"); + endif + + [a, b, c, d, tsam, scaled] = ssdata (sys); + dt = isdt (sys); + + ## default arguments + beta = 1; # ? + tol1 = 0; + tol2 = 0; + ordsel = 1; + nr = 0; + + job = 1; ## ? + + if (dt) # discrete-time + alpha = 1; # ALPHA <= 0 + else # continuous-time + alpha = 0; # 0 <= ALPHA <= 1 + endif + + for k = 1 : 2 : (nargin-1) + prop = lower (varargin{k}); + val = varargin{k+1}; + switch (prop) + case {"order", "n", "nr"} + if (! issample (val, 0) || val != round (val)) + error ("bstmodred: argument %s must be an integer >= 0", varargin{k}); + endif + nr = val; + ordsel = 0; + + case "tol1" + if (! is_real_scalar (val)) + error ("hnamodred: argument %s must be a real scalar", varargin{k}); + endif + tol1 = val; + + case "tol2" + if (! is_real_scalar (val)) + error ("hnamodred: argument %s must be a real scalar", varargin{k}); + endif + tol2 = val; + + case "alpha" + if (! is_real_scalar (val)) + error ("bstmodred: argument %s must be a real scalar", varargin{k}); + endif + if (dt) # discrete-time + if (val < 0 || val > 1) + error ("bstmodred: argument %s must be 0 <= ALPHA <= 1", varargin{k}); + endif + else # continuous-time + if (val > 0) + error ("bstmodred: argument %s must be ALPHA <= 0", varargin{k}); + endif + endif + alpha = val; + + case "beta" + if (! issample (val, 0)) + error ("bstmodred: argument %s must be BETA >= 0", varargin{k}); + endif + beta = val; + + otherwise + error ("hnamodred: invalid property name"); + endswitch + endfor + + ## TODO: handle job + + ## perform model order reduction + [ar, br, cr, dr, nr] = slab09hd (a, b, c, d, dt, scaled, job, nr, ordsel, alpha, beta, \ + tol1, tol2); + + ## assemble reduced order model + sysr = ss (ar, br, cr, dr, tsam); + +endfunction + + +%!shared Mo, Me +%! A = [ -0.04165 0.0000 4.9200 -4.9200 0.0000 0.0000 0.0000 +%! -5.2100 -12.500 0.0000 0.0000 0.0000 0.0000 0.0000 +%! 0.0000 3.3300 -3.3300 0.0000 0.0000 0.0000 0.0000 +%! 0.5450 0.0000 0.0000 0.0000 -0.5450 0.0000 0.0000 +%! 0.0000 0.0000 0.0000 4.9200 -0.04165 0.0000 4.9200 +%! 0.0000 0.0000 0.0000 0.0000 -5.2100 -12.500 0.0000 +%! 0.0000 0.0000 0.0000 0.0000 0.0000 3.3300 -3.3300 ]; +%! +%! B = [ 0.0000 0.0000 +%! 12.500 0.0000 +%! 0.0000 0.0000 +%! 0.0000 0.0000 +%! 0.0000 0.0000 +%! 0.0000 12.500 +%! 0.0000 0.0000 ]; +%! +%! C = [ 1.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 +%! 0.0000 0.0000 0.0000 1.0000 0.0000 0.0000 0.0000 +%! 0.0000 0.0000 0.0000 0.0000 1.0000 0.0000 0.0000 ]; +%! +%! D = [ 0.0000 0.0000 +%! 0.0000 0.0000 +%! 0.0000 0.0000 ]; +%! +%! sys = ss (A, B, C, D, "scaled", true); +%! +%! sysr = bstmodred (sys, "beta", 1.0, "tol1", 0.1, "tol2", 0.0); +%! [Ao, Bo, Co, Do] = ssdata (sysr); +%! +%! Ae = [ 1.2729 0.0000 6.5947 0.0000 -3.4229 +%! 0.0000 0.8169 0.0000 2.4821 0.0000 +%! -2.9889 0.0000 -2.9028 0.0000 -0.3692 +%! 0.0000 -3.3921 0.0000 -3.1126 0.0000 +%! -1.4767 0.0000 -2.0339 0.0000 -0.6107 ]; +%! +%! Be = [ 0.1331 -0.1331 +%! -0.0862 -0.0862 +%! -2.6777 2.6777 +%! -3.5767 -3.5767 +%! -2.3033 2.3033 ]; +%! +%! Ce = [ -0.6907 -0.6882 0.0779 0.0958 -0.0038 +%! 0.0676 0.0000 0.6532 0.0000 -0.7522 +%! 0.6907 -0.6882 -0.0779 0.0958 0.0038 ]; +%! +%! De = [ 0.0000 0.0000 +%! 0.0000 0.0000 +%! 0.0000 0.0000 ]; +%! +%! Mo = [Ao, Bo; Co, Do]; +%! Me = [Ae, Be; Ce, De]; +%! +%!assert (Mo, Me, 1e-4); + Modified: trunk/octave-forge/extra/control-devel/inst/hnamodred.m =================================================================== --- trunk/octave-forge/extra/control-devel/inst/hnamodred.m 2011-10-26 18:59:06 UTC (rev 8867) +++ trunk/octave-forge/extra/control-devel/inst/hnamodred.m 2011-10-26 19:13:06 UTC (rev 8868) @@ -71,7 +71,7 @@ tol1 = 0; tol2 = 0; ordsel = 1; - nr = 0 + nr = 0; if (dt) # discrete-time alpha = 1; # ALPHA <= 0 This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2011-10-27 06:04:38
|
Revision: 8869 http://octave.svn.sourceforge.net/octave/?rev=8869&view=rev Author: paramaniac Date: 2011-10-27 06:04:32 +0000 (Thu, 27 Oct 2011) Log Message: ----------- control-devel: add draft code Modified Paths: -------------- trunk/octave-forge/extra/control-devel/INDEX Added Paths: ----------- trunk/octave-forge/extra/control-devel/src/slab09id.cc Modified: trunk/octave-forge/extra/control-devel/INDEX =================================================================== --- trunk/octave-forge/extra/control-devel/INDEX 2011-10-26 19:13:06 UTC (rev 8868) +++ trunk/octave-forge/extra/control-devel/INDEX 2011-10-27 06:04:32 UTC (rev 8869) @@ -6,9 +6,13 @@ @iddata/set @iddata/size System Identification + fitfrd moesp n4sid Model Reduction - hna + bstmodred + btamodred + hnamodred + spamodred Controller Reduction Added: trunk/octave-forge/extra/control-devel/src/slab09id.cc =================================================================== --- trunk/octave-forge/extra/control-devel/src/slab09id.cc (rev 0) +++ trunk/octave-forge/extra/control-devel/src/slab09id.cc 2011-10-27 06:04:32 UTC (rev 8869) @@ -0,0 +1,425 @@ +/* + +Copyright (C) 2011 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/>. + +Model reduction based on Balance & Truncate (B&T) or +Singular Perturbation Approximation (SPA) method. +Uses SLICOT AB09ID by courtesy of NICONET e.V. +<http://www.slicot.org> + +Author: Lukas Reichlin <luk...@gm...> +Created: October 2011 +Version: 0.1 + +*/ + +#include <octave/oct.h> +#include <f77-fcn.h> +#include "common.cc" + +extern "C" +{ + int F77_FUNC (ab09id, AB09ID) + (char& JOBV, char& JOBW, char& JOBINV, + char& DICO, char& EQUIL, char& ORDSEL, + int& N, int& NV, int& NW, int& M, int& P, + int& NR, + double& ALPHA, + double* A, int& LDA, + double* B, int& LDB, + double* C, int& LDC, + double* D, int& LDD, + double* AV, int& LDAV, + double* BV, int& LDBV, + double* CV, int& LDCV, + double* DV, int& LDDV, + double* AW, int& LDAW, + double* BW, int& LDBW, + double* CW, int& LDCW, + double* DW, int& LDDW, + int& NS, + double* HSV, + double& TOL1, double& TOL2, + int* IWORK, + double* DWORK, int& LDWORK, + int& IWARN, int& INFO); +} + +DEFUN_DLD (slab09id, args, nargout, + "-*- texinfo -*-\n\ +Slicot AB09ID Release 5.0\n\ +No argument checking.\n\ +For internal use only.") +{ + int nargin = args.length (); + octave_value_list retval; + + if (nargin != 22) + { + print_usage (); + } + else + { + // arguments in + char jobv; + char jobw; + char jobinv; + char dico; + char equil; + char ordsel; + + Matrix a = args(0).matrix_value (); + Matrix b = args(1).matrix_value (); + Matrix c = args(2).matrix_value (); + Matrix d = args(3).matrix_value (); + + const int idico = args(4).int_value (); + const int iequil = args(5).int_value (); + int nr = args(6).int_value (); + const int iordsel = args(7).int_value (); + double alpha = args(8).double_value (); + + const int ijobv = args(9).int_value (); + Matrix av = args(10).matrix_value (); + Matrix bv = args(11).matrix_value (); + Matrix cv = args(12).matrix_value (); + Matrix dv = args(13).matrix_value (); + + const int ijobw = args(14).int_value (); + Matrix aw = args(15).matrix_value (); + Matrix bw = args(16).matrix_value (); + Matrix cw = args(17).matrix_value (); + Matrix dw = args(18).matrix_value (); + + const int ijobinv = args(19).int_value (); + double tol1 = args(20).double_value (); + double tol2 = args(21).double_value (); + + switch (ijobv) + { + case 0: + jobv = 'N'; + break; + case 1: + jobv = 'V'; + break; + case 2: + jobv = 'I'; + break; + case 3: + jobv = 'C'; + break; + case 4: + jobv = 'R'; + break; + default: + error ("slab09id: argument jobv invalid"); + } + + switch (ijobw) + { + case 0: + jobw = 'N'; + break; + case 1: + jobw = 'W'; + break; + case 2: + jobw = 'I'; + break; + case 3: + jobw = 'C'; + break; + case 4: + jobw = 'R'; + break; + default: + error ("slab09id: argument jobw invalid"); + } + + switch (ijobinv) + { + case 0: + jobinv = 'N'; + break; + case 1: + jobinv = 'I'; + break; + case 2: + jobinv = 'A'; + break; + default: + error ("slab09id: argument jobinv invalid"); + } + + if (idico == 0) + dico = 'C'; + else + dico = 'D'; + + if (iequil == 0) + equil = 'S'; + else + equil = 'N'; + + if (iordsel == 0) + ordsel = 'F'; + else + ordsel = 'A'; + + int n = a.rows (); // n: number of states + int nv = av.rows (); + int nw = aw.rows (); + int m = b.columns (); // m: number of inputs + int p = c.rows (); // p: number of outputs + + int lda = max (1, n); + int ldb = max (1, n); + int ldc = max (1, p); + int ldd = max (1, p); + + int ldav = max (1, nv); + int ldbv = max (1, nv); + int ldcv = max (1, p); + int lddv = max (1, p); + + int ldaw = max (1, nw); + int ldbw = max (1, nw); + int ldcw = max (1, m); + int lddw = max (1, m); + + // arguments out + int ns; + ColumnVector hsv (n); + + // workspace + int liwork; + int tmpc; + int tmpd; + + if (jobv == 'N') + tmpc = 0; + else + tmpc = max (2*p, nv+p+n+6, 2*nv+p+2); + + if (jobw == 'N') + tmpd = 0; + else + tmpd = max (2*m, nw+m+n+6, 2*nw+m+2); + + if (dico == 'C') + liwork = max (1, m, tmpc, tmpd); + else + liwork = max (1, n, m, tmpc, tmpd); + + int ldwork; + int nvp = nv + p; + int nwm = nw + m; + int ldw1; + int ldw2; + int ldw3 = n*(2*n + max (n, m, p) + 5) + n*(n+1)/2; + int ldw4 = n*(m+p+2) + 2*m*p + min (n, m) + max (3*m+1, min (n, m) + p); + + if (jobv == 'N') + { + ldw1 = 0; + } + else + { + ldw1 = 2*nvp*(nvp+p) + p*p + max (2*nvp*nvp + max (11*nvp+16, p*nvp), + nvp*n + max (nvp*n+n*n, p*n, p*m)); + } + + if (jobw == 'N') + { + ldw2 = 0; + } + else + { + ldw2 = 2*nwm*(nwm+m) + m*m + max (2*nwm*nwm + max (11*nwm+16, m*nwm), + nwm*n + max (nwm*n+n*n, m*n, p*m)); + } + + ldwork = max (ldw1, ldw2, ldw3, ldw4); + + OCTAVE_LOCAL_BUFFER (int, iwork, liwork); + OCTAVE_LOCAL_BUFFER (double, dwork, ldwork); + + // error indicators + int iwarn = 0; + int info = 0; + + + // SLICOT routine AB09ID + F77_XFCN (ab09id, AB09ID, + (jobv, jobw, jobinv, + dico, equil, ordsel, + n, nv, nw, m, p, + nr, + alpha, + a.fortran_vec (), lda, + b.fortran_vec (), ldb, + c.fortran_vec (), ldc, + d.fortran_vec (), ldd, + av.fortran_vec (), ldav, + bv.fortran_vec (), ldbv, + cv.fortran_vec (), ldcv, + dv.fortran_vec (), lddv, + aw.fortran_vec (), ldaw, + bw.fortran_vec (), ldbw, + cw.fortran_vec (), ldcw, + dw.fortran_vec (), lddw, + ns, + hsv.fortran_vec (), + tol1, tol2, + iwork, + dwork, ldwork, + iwarn, info)); + + if (f77_exception_encountered) + error ("hnamodred: exception in SLICOT subroutine AB09ID"); + + if (info != 0) + { + //error ("hsvd: slab09id: AB09ID returned info = %d", info); + if (info < 0) + error ("hnamodred: the %d-th argument had an invalid value", info); + else + { + switch (info) + { + // FIXME: The code below looks nice, but the error message does not + // because there is much white space after each line break + case 1: + error ("hnamodred: 1: the computation of the ordered real Schur form of A\ + failed"); + case 2: + error ("hnamodred: 2: the separation of the ALPHA-stable/unstable\ + diagonal blocks failed because of very close eigenvalues"); + case 3: + error ("hnamodred: 3: the reduction of AV to a real Schur form failed"); + case 4: + error ("hnamodred: 4: the reduction of AW to a real Schur form failed"); + case 5: + error ("hnamodred: 5: the reduction to generalized Schur form of the\ + descriptor pair corresponding to the inverse of V\ + failed"); + case 6: + error ("hnamodred: 6: the reduction to generalized Schur form of the\ + descriptor pair corresponding to the inverse of W\ + failed"); + case 7: + error ("hnamodred: 7: the computation of Hankel singular values failed"); + case 8: + error ("hnamodred: 8: the computation of stable projection in the\ + Hankel-norm approximation algorithm failed"); + case 9: + error ("hnamodred: 9: the order of computed stable projection in the\ + Hankel-norm approximation algorithm differs\ + from the order of Hankel-norm approximation"); + case 10: + error ("hnamodred: 10: the reduction of AV-BV*inv(DV)*CV to a\ + real Schur form failed"); + case 11: + error ("hnamodred: 11: the reduction of AW-BW*inv(DW)*CW to a\ + real Schur form failed"); + case 12: + error ("hnamodred: 12: the solution of the Sylvester equation failed\ + because the poles of V (if JOBV = 'V') or of\ + conj(V) (if JOBV = 'C') are not distinct from\ + the poles of G1 (see METHOD)"); + case 13: + error ("hnamodred: 13: the solution of the Sylvester equation failed\ + because the poles of W (if JOBW = 'W') or of\ + conj(W) (if JOBW = 'C') are not distinct from\ + the poles of G1 (see METHOD)"); + case 14: + error ("hnamodred: 14: the solution of the Sylvester equation failed\ + because the zeros of V (if JOBV = 'I') or of\ + conj(V) (if JOBV = 'R') are not distinct from\ + the poles of G1sr (see METHOD)"); + case 15: + error ("hnamodred: 15: the solution of the Sylvester equation failed\ + because the zeros of W (if JOBW = 'I') or of\ + conj(W) (if JOBW = 'R') are not distinct from\ + the poles of G1sr (see METHOD)"); + case 16: + error ("hnamodred: 16: the solution of the generalized Sylvester system\ + failed because the zeros of V (if JOBV = 'I') or\ + of conj(V) (if JOBV = 'R') are not distinct from\ + the poles of G1sr (see METHOD)"); + case 17: + error ("hnamodred: 17: the solution of the generalized Sylvester system\ + failed because the zeros of W (if JOBW = 'I') or\ + of conj(W) (if JOBW = 'R') are not distinct from\ + the poles of G1sr (see METHOD)"); + case 18: + error ("hnamodred: 18: op(V) is not antistable"); + case 19: + error ("hnamodred: 19: op(W) is not antistable"); + case 20: + error ("hnamodred: 20: V is not invertible"); + case 21: + error ("hnamodred: 21: W is not invertible"); + default: + error ("hnamodred: unknown error, info = %d", info); + } + } + } + + if (iwarn != 0) + { + switch (iwarn) + { + case 1: + warning ("hnamodred: 1: with ORDSEL = 'F', the selected order NR is greater\ + than NSMIN, the sum of the order of the\ + ALPHA-unstable part and the order of a minimal\ + realization of the ALPHA-stable part of the given\ + system. In this case, the resulting NR is set equal\ + to NSMIN."); + break; + case 2: + warning ("hnamodred: 2: with ORDSEL = 'F', the selected order NR is less\ + than the order of the ALPHA-unstable part of the\ + given system. In this case NR is set equal to the\ + order of the ALPHA-unstable part."); + break; + default: + warning ("hnamodred: unknown warning, iwarn = %d", info); + } + } + + // resize + a.resize (nr, nr); + b.resize (nr, m); + c.resize (p, nr); + hsv.resize (ns); + + // return values + retval(0) = a; + retval(1) = b; + retval(2) = c; + retval(3) = d; + retval(4) = octave_value (nr); + // retval(0) = hsv; + // retval(1) = octave_value (ns); + } + + return retval; +} This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2011-10-27 18:42:57
|
Revision: 8887 http://octave.svn.sourceforge.net/octave/?rev=8887&view=rev Author: paramaniac Date: 2011-10-27 18:42:50 +0000 (Thu, 27 Oct 2011) Log Message: ----------- control-devel: update draft code Modified Paths: -------------- trunk/octave-forge/extra/control-devel/devel/makefile_modred.m trunk/octave-forge/extra/control-devel/src/slab09id.cc Modified: trunk/octave-forge/extra/control-devel/devel/makefile_modred.m =================================================================== --- trunk/octave-forge/extra/control-devel/devel/makefile_modred.m 2011-10-27 17:07:19 UTC (rev 8886) +++ trunk/octave-forge/extra/control-devel/devel/makefile_modred.m 2011-10-27 18:42:50 UTC (rev 8887) @@ -3,7 +3,18 @@ srcdir = [develdir, "/../src"]; cd (srcdir); -mkoctfile AB09ID.f TB01PD.f SB08DD.f TB01ID.f TB01KD.f \ +mkoctfile "-Wl,-framework" "-Wl,vecLib" \ + slab09hd.cc \ + AB09HD.f TB01ID.f AB04MD.f TB01KD.f AB09HY.f \ + AB09IX.f MB03UD.f SB02MD.f AB09DD.f TB01LD.f \ + SB03OU.f MA02AD.f MB03QX.f select.f SB03OT.f \ + SB02MR.f SB02MS.f MB03QD.f SB02MU.f SB02MV.f \ + SB02MW.f MB04ND.f MB04OD.f MB03QY.f SB03OR.f \ + SB03OY.f SB04PX.f MB04NY.f MB04OY.f SB03OV.f + +mkoctfile "-Wl,-framework" "-Wl,vecLib" \ + slab09id.cc \ + AB09ID.f TB01PD.f SB08DD.f TB01ID.f TB01KD.f \ AB09IX.f AB09IY.f SB08CD.f MB04ND.f TB01XD.f \ MB04OD.f MB01WD.f MB03UD.f AB07MD.f SB01FY.f \ AB09DD.f TB01LD.f SB03OU.f TB01UD.f MA02AD.f \ @@ -23,14 +34,6 @@ MB04ND.f MB04OD.f SB03OR.f SB03OY.f MB04NY.f \ MB04OY.f SB03OV.f -mkoctfile "-Wl,-framework" "-Wl,vecLib" \ - slab09hd.cc \ - AB09HD.f TB01ID.f AB04MD.f TB01KD.f AB09HY.f \ - AB09IX.f MB03UD.f SB02MD.f AB09DD.f TB01LD.f \ - SB03OU.f MA02AD.f MB03QX.f select.f SB03OT.f \ - SB02MR.f SB02MS.f MB03QD.f SB02MU.f SB02MV.f \ - SB02MW.f MB04ND.f MB04OD.f MB03QY.f SB03OR.f \ - SB03OY.f SB04PX.f MB04NY.f MB04OY.f SB03OV.f system ("rm *.o"); cd (homedir); Modified: trunk/octave-forge/extra/control-devel/src/slab09id.cc =================================================================== --- trunk/octave-forge/extra/control-devel/src/slab09id.cc 2011-10-27 17:07:19 UTC (rev 8886) +++ trunk/octave-forge/extra/control-devel/src/slab09id.cc 2011-10-27 18:42:50 UTC (rev 8887) @@ -81,7 +81,7 @@ char jobc; char jobo; char job; - char weigth + char weight; char equil; char ordsel; @@ -182,7 +182,7 @@ int p = c.rows (); // p: number of outputs int nv = av.rows (); - int pw = cv.rows (); + int pv = cv.rows (); int nw = aw.rows (); int mw = bw.columns (); @@ -233,78 +233,53 @@ else liwrk3 = nw + max (m, mw); - liwork = max (3, liwrk1, liwrk2, liwrk3)); + liwork = max (3, liwrk1, liwrk2, liwrk3); int ldwork; - int - - ldwork = max (lminl, lminr, lrcf, - 2*n*n + max (1, lleft, lright, 2*n*n+5*n, n*max (m, p))); -c where -c lminl = 0, if weight = 'r' or 'n' or nv = 0; otherwise, -c lminl = max(llcf,nv+max(nv,3*p)) if p = pv; -c lminl = max(p,pv)*(2*nv+max(p,pv))+ -c max(llcf,nv+max(nv,3*p,3*pv)) if p <> pv; -c lrcf = 0, and -c lminr = 0, if weight = 'l' or 'n' or nw = 0; otherwise, -c lminr = nw+max(nw,3*m) if m = mw; -c lminr = 2*nw*max(m,mw)+nw+max(nw,3*m,3*mw) if m <> mw; -c llcf = pv*(nv+pv)+pv*nv+max(nv*(nv+5), pv*(pv+2), -c 4*pv, 4*p); -c lrcf = mw*(nw+mw)+max(nw*(nw+5),mw*(mw+2),4*mw,4*m) -c lleft = (n+nv)*(n+nv+max(n+nv,pv)+5) -c if weight = 'l' or 'b' and pv > 0; -c lleft = n*(p+5) if weight = 'r' or 'n' or pv = 0; -c lright = (n+nw)*(n+nw+max(n+nw,mw)+5) -c if weight = 'r' or 'b' and mw > 0; -c lright = n*(m+5) if weight = 'l' or 'n' or mw = 0. + int lminl; + int lrcf; + int lminr; + int llcf; + int lleft; + int lright; - - if (jobv == 'N') - tmpc = 0; - else - tmpc = max (2*p, nv+p+n+6, 2*nv+p+2); - - if (jobw == 'N') - tmpd = 0; - else - tmpd = max (2*m, nw+m+n+6, 2*nw+m+2); - - if (dico == 'C') - liwork = max (1, m, tmpc, tmpd); - else - liwork = max (1, n, m, tmpc, tmpd); - - int ldwork; - int nvp = nv + p; - int nwm = nw + m; - int ldw1; - int ldw2; - int ldw3 = n*(2*n + max (n, m, p) + 5) + n*(n+1)/2; - int ldw4 = n*(m+p+2) + 2*m*p + min (n, m) + max (3*m+1, min (n, m) + p); - - if (jobv == 'N') + if (nw == 0 || weight == 'L' || weight == 'N') { - ldw1 = 0; + lrcf = 0; + lminr = 0; } else { - ldw1 = 2*nvp*(nvp+p) + p*p + max (2*nvp*nvp + max (11*nvp+16, p*nvp), - nvp*n + max (nvp*n+n*n, p*n, p*m)); + lrcf = mw*(nw+mw) + max (nw*(nw+5), mw*(mw+2), 4*mw, 4*m); + if (m == mw) + lminr = nw + max (nw, 3*m); + else + lminr = 2*nw*max (m, mw) + nw + max (nw, 3*m, 3*mw); } - if (jobw == 'N') - { - ldw2 = 0; - } + llcf = pv*(nv+pv) + pv*nv + max (nv*(nv+5), pv*(pv+2), 4*pv, 4*p); + + if (nv == 0 || weight == 'R' || weight == 'N') + lminl = 0; + else if (p == pv) + lminl = max (llcf, nv + max (nv, 3*p)); else - { - ldw2 = 2*nwm*(nwm+m) + m*m + max (2*nwm*nwm + max (11*nwm+16, m*nwm), - nwm*n + max (nwm*n+n*n, m*n, p*m)); - } - - ldwork = max (ldw1, ldw2, ldw3, ldw4); + lminl = max (p, pv) * (2*nv + max (p, pv)) + max (llcf, nv + max (nv, 3*p, 3*pv)); + + if (pv == 0 || weight == 'R' || weight == 'N') + lleft = n*(p+5); + else + lleft = (n+nv) * (n + nv + max (n+nv, pv) + 5); + + if (mw == 0 || weight == 'L' || weight == 'N') + lright = n*(m+5); + else + lright = (n+nw) * (n + nw + max (n+nw, mw) + 5); + + ldwork = max (lminl, lminr, lrcf, + 2*n*n + max (1, lleft, lright, 2*n*n+5*n, n*max (m, p))); + OCTAVE_LOCAL_BUFFER (int, iwork, liwork); OCTAVE_LOCAL_BUFFER (double, dwork, ldwork); This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2011-10-27 21:44:56
|
Revision: 8889 http://octave.svn.sourceforge.net/octave/?rev=8889&view=rev Author: paramaniac Date: 2011-10-27 21:44:50 +0000 (Thu, 27 Oct 2011) Log Message: ----------- control-devel: get testcase for slicot ab09id working Modified Paths: -------------- trunk/octave-forge/extra/control-devel/src/slab09id.cc Added Paths: ----------- trunk/octave-forge/extra/control-devel/devel/test_ab09id.m Added: trunk/octave-forge/extra/control-devel/devel/test_ab09id.m =================================================================== --- trunk/octave-forge/extra/control-devel/devel/test_ab09id.m (rev 0) +++ trunk/octave-forge/extra/control-devel/devel/test_ab09id.m 2011-10-27 21:44:50 UTC (rev 8889) @@ -0,0 +1,82 @@ +% AB09ID EXAMPLE PROGRAM DATA (Continuous system) +% 3 1 1 6 1 0 0 2 0.0 0.0 0.0 0.1E0 0.0 C S S F L S F + +% N, M, P, NV, PV, NW, MW, NR, +% ALPHA, ALPHAC, ALPHAO, TOL1, TOL2, +% DICO, JOBC, JOBO, JOB, WEIGHT, +% EQUIL, ORDSEL + + +a = [ -26.4000, 6.4023, 4.3868; + 32.0000, 0, 0; + 0, 8.0000, 0 ]; + +b = [ 16 + 0 + 0 ]; + +c = [ 9.2994 1.1624 0.1090 ]; + +d = [ 0 ]; + +av = [ -1.0000, 0, 4.0000, -9.2994, -1.1624, -0.1090; + 0, 2.0000, 0, -9.2994, -1.1624, -0.1090; + 0, 0, -3.0000, -9.2994, -1.1624, -0.1090; + 16.0000, 16.0000, 16.0000, -26.4000, 6.4023, 4.3868; + 0, 0, 0, 32.0000, 0, 0; + 0, 0, 0, 0, 8.0000, 0 ]; + +bv = [ 1 + 1 + 1 + 0 + 0 + 0 ]; + +cv = [ 1 1 1 0 0 0 ]; + +dv = [ 0 ]; + +aw = bw = cw = dw = []; + +alpha = alphac = alphao = 0.0; +tol1 = 0.1; +tol2 = 0.0; +dico = 0; +jobc = jobo = 0; +job = 1; +weight = 1; +equil = 0; +ordsel = 0; +nr = 2; + + +[ar, br, cr, dr] = slab09id (a, b, c, d, dico, equil, nr, ordsel, alpha, job, \ + av, bv, cv, dv, \ + aw, bw, cw, dw, \ + weight, jobc, jobo, alphac, alphao, \ + tol1, tol2) + +%{ + AB09ID EXAMPLE PROGRAM RESULTS + + + The order of reduced model = 2 + + The Hankel singular values of weighted ALPHA-stable part are + 3.8253 0.2005 + + The reduced state dynamics matrix Ar is + 9.1900 0.0000 + 0.0000 -34.5297 + + The reduced input/state matrix Br is + 11.9593 + 16.9329 + + The reduced state/output matrix Cr is + 2.8955 6.9152 + + The reduced input/output matrix Dr is + 0.0000 +%} \ No newline at end of file Modified: trunk/octave-forge/extra/control-devel/src/slab09id.cc =================================================================== --- trunk/octave-forge/extra/control-devel/src/slab09id.cc 2011-10-27 19:03:15 UTC (rev 8888) +++ trunk/octave-forge/extra/control-devel/src/slab09id.cc 2011-10-27 21:44:50 UTC (rev 8889) @@ -70,7 +70,7 @@ int nargin = args.length (); octave_value_list retval; - if (nargin != 26) + if (nargin != 25) { print_usage (); } @@ -102,19 +102,19 @@ Matrix cv = args(12).matrix_value (); Matrix dv = args(13).matrix_value (); - Matrix aw = args(15).matrix_value (); - Matrix bw = args(16).matrix_value (); - Matrix cw = args(17).matrix_value (); - Matrix dw = args(18).matrix_value (); + Matrix aw = args(14).matrix_value (); + Matrix bw = args(15).matrix_value (); + Matrix cw = args(16).matrix_value (); + Matrix dw = args(17).matrix_value (); - const int iweight = args(19).int_value (); - const int ijobc = args(20).int_value (); - double alphac = args(21).double_value (); - const int ijobo = args(22).int_value (); - double alphao = args(23).double_value (); + const int iweight = args(18).int_value (); + const int ijobc = args(19).int_value (); + double alphac = args(20).double_value (); + const int ijobo = args(21).int_value (); + double alphao = args(22).double_value (); - double tol1 = args(24).double_value (); - double tol2 = args(25).double_value (); + double tol1 = args(23).double_value (); + double tol2 = args(24).double_value (); if (idico == 0) dico = 'C'; This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2011-11-10 17:31:03
|
Revision: 9046 http://octave.svn.sourceforge.net/octave/?rev=9046&view=rev Author: paramaniac Date: 2011-11-10 17:30:57 +0000 (Thu, 10 Nov 2011) Log Message: ----------- control-devel: update makefile Modified Paths: -------------- trunk/octave-forge/extra/control-devel/inst/__check_weight__.m trunk/octave-forge/extra/control-devel/src/Makefile Modified: trunk/octave-forge/extra/control-devel/inst/__check_weight__.m =================================================================== --- trunk/octave-forge/extra/control-devel/inst/__check_weight__.m 2011-11-10 07:14:45 UTC (rev 9045) +++ trunk/octave-forge/extra/control-devel/inst/__check_weight__.m 2011-11-10 17:30:57 UTC (rev 9046) @@ -33,5 +33,7 @@ [a, b, c, d] = ssdata (sys); job = 1; + + ## TODO: check system size endfunction Modified: trunk/octave-forge/extra/control-devel/src/Makefile =================================================================== --- trunk/octave-forge/extra/control-devel/src/Makefile 2011-11-10 07:14:45 UTC (rev 9045) +++ trunk/octave-forge/extra/control-devel/src/Makefile 2011-11-10 17:30:57 UTC (rev 9046) @@ -1,10 +1,10 @@ -all: slab09hd.oct slab09jd.oct slsb10yd.oct +all: slab09hd.oct slab09id.oct slab09jd.oct slsb10yd.oct # TODO: leading and trailing underscores for sl* functions # (__sl*__.oct) would be nice, but this can be an issue # for fortran compilers. -# model reduction based on balanced stochastic truncation method +# balanced stochastic truncation model reduction slab09hd.oct: slab09hd.cc mkoctfile slab09hd.cc \ AB09HD.f TB01ID.f AB04MD.f TB01KD.f AB09HY.f \ @@ -14,7 +14,19 @@ SB02MW.f MB04ND.f MB04OD.f MB03QY.f SB03OR.f \ SB03OY.f SB04PX.f MB04NY.f MB04OY.f SB03OV.f -# frequency-weighted Hankel norm approximation with invertible weights +# balanced truncation & singular perturbation approximation model reduction +slab09id.oct: slab09id.cc + mkoctfile slab09id.cc \ + AB09ID.f TB01PD.f SB08DD.f TB01ID.f TB01KD.f \ + AB09IX.f AB09IY.f SB08CD.f MB04ND.f TB01XD.f \ + MB04OD.f MB01WD.f MB03UD.f AB07MD.f SB01FY.f \ + AB09DD.f TB01LD.f SB03OU.f TB01UD.f MA02AD.f \ + MA02BD.f MB03OY.f MB03QX.f MB01PD.f select.f \ + MB01YD.f MB04NY.f MB01ZD.f SB03OT.f MB04OX.f \ + MB04OY.f MB03QD.f SB03OY.f MB03QY.f MB01QD.f \ + SB03OR.f SB03OV.f SB04PX.f + +# hankel-norm approximation model reduction slab09jd.oct: slab09jd.cc mkoctfile slab09jd.cc \ AB09JD.f TB01ID.f TB01KD.f AB07ND.f AB09JV.f \ This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2011-11-10 20:31:53
|
Revision: 9054 http://octave.svn.sourceforge.net/octave/?rev=9054&view=rev Author: paramaniac Date: 2011-11-10 20:31:47 +0000 (Thu, 10 Nov 2011) Log Message: ----------- control-devel: touch up hnamodred Modified Paths: -------------- trunk/octave-forge/extra/control-devel/inst/hnamodred.m trunk/octave-forge/extra/control-devel/src/slab09jd.cc Modified: trunk/octave-forge/extra/control-devel/inst/hnamodred.m =================================================================== --- trunk/octave-forge/extra/control-devel/inst/hnamodred.m 2011-11-10 20:28:06 UTC (rev 9053) +++ trunk/octave-forge/extra/control-devel/inst/hnamodred.m 2011-11-10 20:31:47 UTC (rev 9054) @@ -16,7 +16,7 @@ ## along with LTI Syncope. If not, see <http://www.gnu.org/licenses/>. ## -*- texinfo -*- -## @deftypefn{Function File} {[@var{sysr}, @var{nr}] =} hnamodred (@var{sys}, @dots{}) +## @deftypefn{Function File} {[@var{sysr}, @var{info}] =} hnamodred (@var{sys}, @dots{}) ## Model order reduction by frequency weighted optimal Hankel-norm approximation method. ## ## @strong{Inputs} @@ -45,7 +45,7 @@ ## Created: October 2011 ## Version: 0.1 -function [sysr, nr] = hnamodred (sys, varargin) +function [sysr, info] = hnamodred (sys, varargin) if (nargin == 0) print_usage (); @@ -87,6 +87,30 @@ case {"right", "w"} [aw, bw, cw, dw, jobw] = __modred_check_weight__ (val, dt, m, m); + case {"left-inv", "inv-v"} + [av, bv, cv, dv] = __modred_check_weight__ (val, dt, p, p); + jobv = 2; + + case {"right-inv", "inv-w"} + [aw, bw, cw, dw] = __modred_check_weight__ (val, dt, m, m); + jobv = 2 + + case {"left-conj", "conj-v"} + [av, bv, cv, dv] = __modred_check_weight__ (val, dt, p, p); + jobv = 3; + + case {"right-conj", "conj-w"} + [aw, bw, cw, dw] = __modred_check_weight__ (val, dt, m, m); + jobv = 3 + + case {"left-conj-inv", "conj-inv-v"} + [av, bv, cv, dv] = __modred_check_weight__ (val, dt, p, p); + jobv = 4; + + case {"right-conj-inv", "conj-inv-w"} + [aw, bw, cw, dw] = __modred_check_weight__ (val, dt, m, m); + jobv = 4 + case {"order", "n", "nr"} [nr, ordsel] = __modred_check_order__ (val); @@ -99,6 +123,18 @@ case "alpha" alpha = __modred_check_alpha__ (val, dt); + case {"approach", "jobinv"} + switch (tolower (val(1))) + case {"d", "n"} # "descriptor" + jobinv = 0; + case {"s", "i"} # "standard" + jobinv = 1; + case "a" # {"auto", "automatic"} + jobinv = 2; + otherwise + error ("hnamodred: invalid computational approach"); + endswitch + otherwise error ("hnamodred: invalid property name"); endswitch @@ -107,14 +143,17 @@ ## TODO: handle jobv, jobw, (jobinv) ## perform model order reduction - [ar, br, cr, dr, nr] = slab09jd (a, b, c, d, dt, scaled, nr, ordsel, alpha, \ - jobv, av, bv, cv, dv, \ - jobw, aw, bw, cw, dw, \ - jobinv, tol1, tol2); + [ar, br, cr, dr, nr, hsv, ns] = slab09jd (a, b, c, d, dt, scaled, nr, ordsel, alpha, \ + jobv, av, bv, cv, dv, \ + jobw, aw, bw, cw, dw, \ + jobinv, tol1, tol2); ## assemble reduced order model sysr = ss (ar, br, cr, dr, tsam); + ## assemble info struct + info = struct ("nr", nr, "ns", ns, "hsv", hsv); + endfunction Modified: trunk/octave-forge/extra/control-devel/src/slab09jd.cc =================================================================== --- trunk/octave-forge/extra/control-devel/src/slab09jd.cc 2011-11-10 20:28:06 UTC (rev 9053) +++ trunk/octave-forge/extra/control-devel/src/slab09jd.cc 2011-11-10 20:31:47 UTC (rev 9054) @@ -414,8 +414,8 @@ retval(2) = c; retval(3) = d; retval(4) = octave_value (nr); - // retval(0) = hsv; - // retval(1) = octave_value (ns); + retval(5) = hsv; + retval(6) = octave_value (ns); } return retval; This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2011-11-10 20:38:48
|
Revision: 9055 http://octave.svn.sourceforge.net/octave/?rev=9055&view=rev Author: paramaniac Date: 2011-11-10 20:38:42 +0000 (Thu, 10 Nov 2011) Log Message: ----------- control-devel: return info struct for all model reduction commands Modified Paths: -------------- trunk/octave-forge/extra/control-devel/inst/__modred_ab09id__.m trunk/octave-forge/extra/control-devel/inst/bstmodred.m trunk/octave-forge/extra/control-devel/inst/btamodred.m trunk/octave-forge/extra/control-devel/inst/spamodred.m trunk/octave-forge/extra/control-devel/src/slab09hd.cc trunk/octave-forge/extra/control-devel/src/slab09id.cc Modified: trunk/octave-forge/extra/control-devel/inst/__modred_ab09id__.m =================================================================== --- trunk/octave-forge/extra/control-devel/inst/__modred_ab09id__.m 2011-11-10 20:31:47 UTC (rev 9054) +++ trunk/octave-forge/extra/control-devel/inst/__modred_ab09id__.m 2011-11-10 20:38:42 UTC (rev 9055) @@ -16,7 +16,7 @@ ## along with LTI Syncope. If not, see <http://www.gnu.org/licenses/>. ## -*- texinfo -*- -## @deftypefn{Function File} {[@var{sysr}, @var{nr}] =} __modred_ab09id__ (@var{method}, @dots{}) +## @deftypefn{Function File} {[@var{sysr}, @var{info}] =} __modred_ab09id__ (@var{method}, @dots{}) ## Backend for btamodred and spamodred. ## @end deftypefn @@ -24,7 +24,7 @@ ## Created: November 2011 ## Version: 0.1 -function [sysr, nr] = __modred_ab09id__ (method, varargin) +function [sysr, info] = __modred_ab09id__ (method, varargin) if (nargin < 2) print_usage (); @@ -110,15 +110,18 @@ ## TODO: handle job ## perform model order reduction - [ar, br, cr, dr] = slab09id (a, b, c, d, dico, equil, nr, ordsel, alpha, job, \ - av, bv, cv, dv, \ - aw, bw, cw, dw, \ - weight, jobc, jobo, alphac, alphao, \ - tol1, tol2); + [ar, br, cr, dr, nr, hsv, ns] = slab09id (a, b, c, d, dico, equil, nr, ordsel, alpha, job, \ + av, bv, cv, dv, \ + aw, bw, cw, dw, \ + weight, jobc, jobo, alphac, alphao, \ + tol1, tol2); ## assemble reduced order model sysr = ss (ar, br, cr, dr, tsam); + ## assemble info struct + info = struct ("nr", nr, "ns", ns, "hsv", hsv); + endfunction Modified: trunk/octave-forge/extra/control-devel/inst/bstmodred.m =================================================================== --- trunk/octave-forge/extra/control-devel/inst/bstmodred.m 2011-11-10 20:31:47 UTC (rev 9054) +++ trunk/octave-forge/extra/control-devel/inst/bstmodred.m 2011-11-10 20:38:42 UTC (rev 9055) @@ -103,12 +103,15 @@ ## TODO: handle job ## perform model order reduction - [ar, br, cr, dr, nr] = slab09hd (a, b, c, d, dt, scaled, job, nr, ordsel, alpha, beta, \ - tol1, tol2); + [ar, br, cr, dr, nr, hsv, ns] = slab09hd (a, b, c, d, dt, scaled, job, nr, ordsel, alpha, beta, \ + tol1, tol2); ## assemble reduced order model sysr = ss (ar, br, cr, dr, tsam); + ## assemble info struct + info = struct ("nr", nr, "ns", ns, "hsv", hsv); + endfunction Modified: trunk/octave-forge/extra/control-devel/inst/btamodred.m =================================================================== --- trunk/octave-forge/extra/control-devel/inst/btamodred.m 2011-11-10 20:31:47 UTC (rev 9054) +++ trunk/octave-forge/extra/control-devel/inst/btamodred.m 2011-11-10 20:38:42 UTC (rev 9055) @@ -45,9 +45,9 @@ ## Created: November 2011 ## Version: 0.1 -function [sysr, nr] = btamodred (varargin) +function [sysr, info] = btamodred (varargin) - [sysr, nr] = __modred_ab09id__ ("bta", varargin{:}); + [sysr, info] = __modred_ab09id__ ("bta", varargin{:}); endfunction Modified: trunk/octave-forge/extra/control-devel/inst/spamodred.m =================================================================== --- trunk/octave-forge/extra/control-devel/inst/spamodred.m 2011-11-10 20:31:47 UTC (rev 9054) +++ trunk/octave-forge/extra/control-devel/inst/spamodred.m 2011-11-10 20:38:42 UTC (rev 9055) @@ -45,9 +45,9 @@ ## Created: November 2011 ## Version: 0.1 -function [sysr, nr] = spamodred (varargin) +function [sysr, info] = spamodred (varargin) - [sysr, nr] = __modred_ab09id__ ("spa", varargin{:}); + [sysr, info] = __modred_ab09id__ ("spa", varargin{:}); endfunction Modified: trunk/octave-forge/extra/control-devel/src/slab09hd.cc =================================================================== --- trunk/octave-forge/extra/control-devel/src/slab09hd.cc 2011-11-10 20:31:47 UTC (rev 9054) +++ trunk/octave-forge/extra/control-devel/src/slab09hd.cc 2011-11-10 20:38:42 UTC (rev 9055) @@ -261,8 +261,8 @@ retval(2) = c; retval(3) = d; retval(4) = octave_value (nr); - // retval(0) = hsv; - // retval(1) = octave_value (ns); + retval(5) = hsv; + retval(6) = octave_value (ns); } return retval; Modified: trunk/octave-forge/extra/control-devel/src/slab09id.cc =================================================================== --- trunk/octave-forge/extra/control-devel/src/slab09id.cc 2011-11-10 20:31:47 UTC (rev 9054) +++ trunk/octave-forge/extra/control-devel/src/slab09id.cc 2011-11-10 20:38:42 UTC (rev 9055) @@ -414,8 +414,8 @@ retval(2) = c; retval(3) = d; retval(4) = octave_value (nr); - // retval(0) = hsv; - // retval(1) = octave_value (ns); + retval(5) = hsv; + retval(6) = octave_value (ns); } return retval; This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2011-11-18 16:53:50
|
Revision: 9134 http://octave.svn.sourceforge.net/octave/?rev=9134&view=rev Author: paramaniac Date: 2011-11-18 16:53:44 +0000 (Fri, 18 Nov 2011) Log Message: ----------- control-devel: add draft code Modified Paths: -------------- trunk/octave-forge/extra/control-devel/inst/bstmodred.m Added Paths: ----------- trunk/octave-forge/extra/control-devel/devel/opt2cell.m trunk/octave-forge/extra/control-devel/devel/options.m Added: trunk/octave-forge/extra/control-devel/devel/opt2cell.m =================================================================== --- trunk/octave-forge/extra/control-devel/devel/opt2cell.m (rev 0) +++ trunk/octave-forge/extra/control-devel/devel/opt2cell.m 2011-11-18 16:53:44 UTC (rev 9134) @@ -0,0 +1,12 @@ +function c = opt2cell (opt) + + if (! isstruct (opt)) + error ("opt2cell: argument must be a struct"); + endif + + key = fieldnames (opt); + val = struct2cell (opt); + + c = [key.'; val.'](:); # reshape to {key1; val1; key2; val2; ...} + +endfunction \ No newline at end of file Added: trunk/octave-forge/extra/control-devel/devel/options.m =================================================================== --- trunk/octave-forge/extra/control-devel/devel/options.m (rev 0) +++ trunk/octave-forge/extra/control-devel/devel/options.m 2011-11-18 16:53:44 UTC (rev 9134) @@ -0,0 +1,8 @@ +function opt = options (varargin) + + key = reshape (varargin(1:2:end-1), [], 1); + val = reshape (varargin(2:2:end), [], 1); + + opt = cell2struct (val, key, 1) + +endfunction \ No newline at end of file Modified: trunk/octave-forge/extra/control-devel/inst/bstmodred.m =================================================================== --- trunk/octave-forge/extra/control-devel/inst/bstmodred.m 2011-11-18 14:23:19 UTC (rev 9133) +++ trunk/octave-forge/extra/control-devel/inst/bstmodred.m 2011-11-18 16:53:44 UTC (rev 9134) @@ -16,8 +16,12 @@ ## along with LTI Syncope. If not, see <http://www.gnu.org/licenses/>. ## -*- texinfo -*- -## @deftypefn{Function File} {[@var{sysr}, @var{nr}] =} hnamodred (@var{sys}, @dots{}) -## Model order reduction by frequency weighted optimal Hankel-norm approximation method. +## @deftypefn{Function File} {[@var{sysr}, @var{nr}] =} bstmodred (@var{sys}, @dots{}) +## Model order reduction by Balanced Stochastic Truncation method. +## Uses the stochastic balancing approach in conjunction with the square-root or +## the balancing-free square-root Balance & Truncate (B&T) +## or Singular Perturbation Approximation (SPA) model reduction +## methods for the ALPHA-stable part of the system. ## ## @strong{Inputs} ## @table @var This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2011-11-19 18:01:34
|
Revision: 9139 http://octave.svn.sourceforge.net/octave/?rev=9139&view=rev Author: paramaniac Date: 2011-11-19 18:01:27 +0000 (Sat, 19 Nov 2011) Log Message: ----------- control-devel: allow struct to pass settings Modified Paths: -------------- trunk/octave-forge/extra/control-devel/inst/__modred_ab09id__.m trunk/octave-forge/extra/control-devel/inst/bstmodred.m trunk/octave-forge/extra/control-devel/inst/hnamodred.m Added Paths: ----------- trunk/octave-forge/extra/control-devel/inst/__opt2cell__.m trunk/octave-forge/extra/control-devel/inst/options.m Removed Paths: ------------- trunk/octave-forge/extra/control-devel/devel/opt2cell.m trunk/octave-forge/extra/control-devel/devel/options.m Deleted: trunk/octave-forge/extra/control-devel/devel/opt2cell.m =================================================================== --- trunk/octave-forge/extra/control-devel/devel/opt2cell.m 2011-11-19 15:43:51 UTC (rev 9138) +++ trunk/octave-forge/extra/control-devel/devel/opt2cell.m 2011-11-19 18:01:27 UTC (rev 9139) @@ -1,12 +0,0 @@ -function c = opt2cell (opt) - - if (! isstruct (opt)) - error ("opt2cell: argument must be a struct"); - endif - - key = fieldnames (opt); - val = struct2cell (opt); - - c = [key.'; val.'](:); # reshape to {key1; val1; key2; val2; ...} - -endfunction \ No newline at end of file Deleted: trunk/octave-forge/extra/control-devel/devel/options.m =================================================================== --- trunk/octave-forge/extra/control-devel/devel/options.m 2011-11-19 15:43:51 UTC (rev 9138) +++ trunk/octave-forge/extra/control-devel/devel/options.m 2011-11-19 18:01:27 UTC (rev 9139) @@ -1,16 +0,0 @@ -function opt = options (varargin) - - if (nargin == 0) - print_usage (); - endif - - if (rem (nargin, 2)) - error ("options: properties and values must come in pairs"); - endif - - key = reshape (varargin(1:2:end-1), [], 1); - val = reshape (varargin(2:2:end), [], 1); - - opt = cell2struct (val, key, 1) - -endfunction \ No newline at end of file Modified: trunk/octave-forge/extra/control-devel/inst/__modred_ab09id__.m =================================================================== --- trunk/octave-forge/extra/control-devel/inst/__modred_ab09id__.m 2011-11-19 15:43:51 UTC (rev 9138) +++ trunk/octave-forge/extra/control-devel/inst/__modred_ab09id__.m 2011-11-19 18:01:27 UTC (rev 9139) @@ -42,7 +42,9 @@ error ("%smodred: first argument must be an LTI system", method); endif - if (rem (npv, 2)) + if (npv == 1) + varargin = __opt2cell__ (varargin{1}); + elseif (rem (npv, 2)) error ("%smodred: properties and values must come in pairs", method); endif @@ -103,7 +105,7 @@ ## TODO: alphac, alphao, jobc, jobo otherwise - error ("hnamodred: invalid property name"); + warning ("modred: invalid property name ""%s"" ignored", prop); endswitch endfor Copied: trunk/octave-forge/extra/control-devel/inst/__opt2cell__.m (from rev 9138, trunk/octave-forge/extra/control-devel/devel/opt2cell.m) =================================================================== --- trunk/octave-forge/extra/control-devel/inst/__opt2cell__.m (rev 0) +++ trunk/octave-forge/extra/control-devel/inst/__opt2cell__.m 2011-11-19 18:01:27 UTC (rev 9139) @@ -0,0 +1,12 @@ +function c = __opt2cell__ (opt) + + if (! isstruct (opt)) + error ("opt2cell: argument must be a struct"); + endif + + key = fieldnames (opt); + val = struct2cell (opt); + + c = [key.'; val.'](:); # reshape to {key1; val1; key2; val2; ...} + +endfunction \ No newline at end of file Modified: trunk/octave-forge/extra/control-devel/inst/bstmodred.m =================================================================== --- trunk/octave-forge/extra/control-devel/inst/bstmodred.m 2011-11-19 15:43:51 UTC (rev 9138) +++ trunk/octave-forge/extra/control-devel/inst/bstmodred.m 2011-11-19 18:01:27 UTC (rev 9139) @@ -59,7 +59,9 @@ error ("bstmodred: first argument must be an LTI system"); endif - if (rem (nargin-1, 2)) + if (nargin == 2) + varargin = __opt2cell__ (varargin{1}); + elseif (rem (nargin-1, 2)) error ("bstmodred: properties and values must come in pairs"); endif @@ -113,12 +115,10 @@ endswitch otherwise - error ("hnamodred: invalid property name"); + warning ("bstmodred: invalid property name ""%s"" ignored", prop); endswitch endfor - ## TODO: handle job - ## perform model order reduction [ar, br, cr, dr, nr, hsv, ns] = slab09hd (a, b, c, d, dt, scaled, job, nr, ordsel, alpha, beta, \ tol1, tol2); Modified: trunk/octave-forge/extra/control-devel/inst/hnamodred.m =================================================================== --- trunk/octave-forge/extra/control-devel/inst/hnamodred.m 2011-11-19 15:43:51 UTC (rev 9138) +++ trunk/octave-forge/extra/control-devel/inst/hnamodred.m 2011-11-19 18:01:27 UTC (rev 9139) @@ -55,7 +55,9 @@ error ("hnamodred: first argument must be an LTI system"); endif - if (rem (nargin-1, 2)) + if (nargin == 2) + varargin = __opt2cell__ (varargin{1}); + elseif (rem (nargin-1, 2)) error ("hnamodred: properties and values must come in pairs"); endif @@ -136,12 +138,11 @@ endswitch otherwise - error ("hnamodred: invalid property name"); + warning ("hnamodred: invalid property name ""%s"" ignored", prop); endswitch endfor + - ## TODO: handle jobv, jobw, (jobinv) - ## perform model order reduction [ar, br, cr, dr, nr, hsv, ns] = slab09jd (a, b, c, d, dt, scaled, nr, ordsel, alpha, \ jobv, av, bv, cv, dv, \ Copied: trunk/octave-forge/extra/control-devel/inst/options.m (from rev 9138, trunk/octave-forge/extra/control-devel/devel/options.m) =================================================================== --- trunk/octave-forge/extra/control-devel/inst/options.m (rev 0) +++ trunk/octave-forge/extra/control-devel/inst/options.m 2011-11-19 18:01:27 UTC (rev 9139) @@ -0,0 +1,16 @@ +function opt = options (varargin) + + if (nargin == 0) + print_usage (); + endif + + if (rem (nargin, 2)) + error ("options: properties and values must come in pairs"); + endif + + key = reshape (varargin(1:2:end-1), [], 1); + val = reshape (varargin(2:2:end), [], 1); + + opt = cell2struct (val, key, 1) + +endfunction \ No newline at end of file This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2011-11-25 16:52:02
|
Revision: 9183 http://octave.svn.sourceforge.net/octave/?rev=9183&view=rev Author: paramaniac Date: 2011-11-25 16:51:56 +0000 (Fri, 25 Nov 2011) Log Message: ----------- control-devel: move file Added Paths: ----------- trunk/octave-forge/extra/control-devel/inst/test_devel.m Removed Paths: ------------- trunk/octave-forge/extra/control-devel/devel/test_devel.m Deleted: trunk/octave-forge/extra/control-devel/devel/test_devel.m =================================================================== --- trunk/octave-forge/extra/control-devel/devel/test_devel.m 2011-11-25 16:46:09 UTC (rev 9182) +++ trunk/octave-forge/extra/control-devel/devel/test_devel.m 2011-11-25 16:51:56 UTC (rev 9183) @@ -1,7 +0,0 @@ -# identification -test fitfrd - -# model order reduction -test bstmodred -test btamodred -test hnamodred Copied: trunk/octave-forge/extra/control-devel/inst/test_devel.m (from rev 9181, trunk/octave-forge/extra/control-devel/devel/test_devel.m) =================================================================== --- trunk/octave-forge/extra/control-devel/inst/test_devel.m (rev 0) +++ trunk/octave-forge/extra/control-devel/inst/test_devel.m 2011-11-25 16:51:56 UTC (rev 9183) @@ -0,0 +1,7 @@ +# identification +test fitfrd + +# model order reduction +test bstmodred +test btamodred +test hnamodred This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2011-11-29 14:00:29
|
Revision: 9220 http://octave.svn.sourceforge.net/octave/?rev=9220&view=rev Author: paramaniac Date: 2011-11-29 14:00:17 +0000 (Tue, 29 Nov 2011) Log Message: ----------- control-devel: update draft code Modified Paths: -------------- trunk/octave-forge/extra/control-devel/devel/makefile_conred.m trunk/octave-forge/extra/control-devel/inst/bstmodred.m trunk/octave-forge/extra/control-devel/src/slsb16ad.cc Modified: trunk/octave-forge/extra/control-devel/devel/makefile_conred.m =================================================================== --- trunk/octave-forge/extra/control-devel/devel/makefile_conred.m 2011-11-29 09:53:30 UTC (rev 9219) +++ trunk/octave-forge/extra/control-devel/devel/makefile_conred.m 2011-11-29 14:00:17 UTC (rev 9220) @@ -12,7 +12,9 @@ srcdir = [develdir, "/../src"]; cd (srcdir); -mkoctfile SB16AD.f TB01ID.f SB16AY.f TB01KD.f AB09IX.f \ +mkoctfile "-Wl,-framework" "-Wl,vecLib" \ + slsb16ad.cc \ + SB16AD.f TB01ID.f SB16AY.f TB01KD.f AB09IX.f \ MB04OD.f MB01WD.f SB03OD.f MB03UD.f AB05PD.f \ AB09DD.f AB07ND.f TB01LD.f AB05QD.f SB03OU.f \ MA02AD.f MB03QX.f select.f MB01YD.f MB01ZD.f \ Modified: trunk/octave-forge/extra/control-devel/inst/bstmodred.m =================================================================== --- trunk/octave-forge/extra/control-devel/inst/bstmodred.m 2011-11-29 09:53:30 UTC (rev 9219) +++ trunk/octave-forge/extra/control-devel/inst/bstmodred.m 2011-11-29 14:00:17 UTC (rev 9220) @@ -97,9 +97,10 @@ ## @item ## Guaranteed a priori error bound ## @iftex -## @math{|| G^{-1} (G-G_r) ||_{\\infty} <= } +## @tex +## $$ || G^{-1} (G-G_r) ||_{\\infty} \\leq 2 \\sum_{j=r+1}^{n} \\frac{1+\\sigma_j}{1-\\sigma_j} - 1 $$ +## @end tex ## @end iftex -## ## @end itemize ## ## Modified: trunk/octave-forge/extra/control-devel/src/slsb16ad.cc =================================================================== --- trunk/octave-forge/extra/control-devel/src/slsb16ad.cc 2011-11-29 09:53:30 UTC (rev 9219) +++ trunk/octave-forge/extra/control-devel/src/slsb16ad.cc 2011-11-29 14:00:17 UTC (rev 9220) @@ -87,7 +87,7 @@ const int idico = args(4).int_value (); const int iequil = args(5).int_value (); - int nr = args(6).int_value (); + int ncr = args(6).int_value (); const int iordsel = args(7).int_value (); double alpha = args(8).double_value (); const int ijobmr = args(9).int_value (); @@ -210,50 +210,26 @@ liwork = max (1, liwrk1, liwrk2); int ldwork; - int lminl; - int lrcf; - int lminr; - int llcf; - int lleft; - int lright; + int lfreq; + int lsqred; - if (nw == 0 || weight == 'L' || weight == 'N') + if (weight == 'N') { - lrcf = 0; - lminr = 0; + if (equil == 'N') // if WEIGHT = 'N' and EQUIL = 'N' + lfreq = nc*(max (m, p) + 5); + else // if WEIGHT = 'N' and EQUIL = 'S' + lfreq = max (n, nc*(max (m, p) + 5)); + } - else + else // if WEIGHT = 'I' or 'O' or 'P' { - lrcf = mw*(nw+mw) + max (nw*(nw+5), mw*(mw+2), 4*mw, 4*m); - if (m == mw) - lminr = nw + max (nw, 3*m); - else - lminr = 2*nw*max (m, mw) + nw + max (nw, 3*m, 3*mw); + lfreq = (n+nc)*(n+nc+2*m+2*p) + + max ((n+nc)*(n+nc+max(n+nc,m,p)+7), (m+p)*(m+p+4)); } - llcf = pv*(nv+pv) + pv*nv + max (nv*(nv+5), pv*(pv+2), 4*pv, 4*p); + lsqred = max (1, 2*nc*nc+5*nc); + ldwork = 2*nc*nc + max (1, lfreq, lsqred); - if (nv == 0 || weight == 'R' || weight == 'N') - lminl = 0; - else if (p == pv) - lminl = max (llcf, nv + max (nv, 3*p)); - else - lminl = max (p, pv) * (2*nv + max (p, pv)) + max (llcf, nv + max (nv, 3*p, 3*pv)); - - - if (pv == 0 || weight == 'R' || weight == 'N') - lleft = n*(p+5); - else - lleft = (n+nv) * (n + nv + max (n+nv, pv) + 5); - - if (mw == 0 || weight == 'L' || weight == 'N') - lright = n*(m+5); - else - lright = (n+nw) * (n + nw + max (n+nw, mw) + 5); - - ldwork = max (lminl, lminr, lrcf, - 2*n*n + max (1, lleft, lright, 2*n*n+5*n, n*max (m, p))); - OCTAVE_LOCAL_BUFFER (int, iwork, liwork); OCTAVE_LOCAL_BUFFER (double, dwork, ldwork); @@ -353,19 +329,19 @@ } // resize - a.resize (nr, nr); - b.resize (nr, m); - c.resize (p, nr); - hsv.resize (ns); + ac.resize (ncr, ncr); + bc.resize (ncr, p); // p: number of plant outputs + cc.resize (m, ncr); // m: number of plant inputs + hsvc.resize (ncs); // return values - retval(0) = a; - retval(1) = b; - retval(2) = c; - retval(3) = d; - retval(4) = octave_value (nr); - retval(5) = hsv; - retval(6) = octave_value (ns); + retval(0) = ac; + retval(1) = bc; + retval(2) = cc; + retval(3) = dc; + retval(4) = octave_value (ncr); + retval(5) = hsvc; + retval(6) = octave_value (ncs); } return retval; This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2011-11-29 14:15:43
|
Revision: 9221 http://octave.svn.sourceforge.net/octave/?rev=9221&view=rev Author: paramaniac Date: 2011-11-29 14:15:32 +0000 (Tue, 29 Nov 2011) Log Message: ----------- control-devel: get slicot sb16ad working Modified Paths: -------------- trunk/octave-forge/extra/control-devel/src/slsb16ad.cc Added Paths: ----------- trunk/octave-forge/extra/control-devel/devel/test_sb16ad.m Added: trunk/octave-forge/extra/control-devel/devel/test_sb16ad.m =================================================================== --- trunk/octave-forge/extra/control-devel/devel/test_sb16ad.m (rev 0) +++ trunk/octave-forge/extra/control-devel/devel/test_sb16ad.m 2011-11-29 14:15:32 UTC (rev 9221) @@ -0,0 +1,88 @@ +%{ + SB16AD EXAMPLE PROGRAM DATA (Continuous system) + 3 1 1 3 2 0.0 0.1E0 0.0 C S S F I N F + + N, M, P, NC, NCR, ALPHA, TOL1, TOL2, DICO, + JOBC, JOBO, JOBMR, WEIGHT, EQUIL, ORDSEL + + -1. 0. 4. + 0. 2. 0. + 0. 0. -3. + 1. + 1. + 1. + 1. 1. 1. + 0. + -26.4000 6.4023 4.3868 + 32.0000 0 0 + 0 8.0000 0 + -16 + 0 + 0 + 9.2994 1.1624 0.1090 + 0 +%} + +a = [ -1. 0. 4. + 0. 2. 0. + 0. 0. -3. ]; + +b = [ 1. + 1. + 1. ]; + +c = [ 1. 1. 1. ]; + +d = [ 0. ]; + +ac = [ -26.4000, 6.4023, 4.3868; + 32.0000, 0, 0; + 0, 8.0000, 0 ]; + +bc = [ -16 + 0 + 0 ]; + +cc = [ 9.2994 1.1624 0.1090 ]; + +dc = [ 0 ]; + + +alpha = 0.0; +tol1 = 0.1; +tol2 = 0.0; +dico = 0; +jobc = jobo = 0; +jobmr = 1; +weight = 2; +equil = 1; +ordsel = 0; +ncr = 2; + + +[ar, br, cr, dr] = slsb16ad (a, b, c, d, dico, equil, ncr, ordsel, alpha, jobmr, \ + ac, bc, cc, dc, weight, jobc, jobo, tol1, tol2) + +%{ + SB16AD EXAMPLE PROGRAM RESULTS + + + The order of reduced controller = 2 + + The Hankel singular values of weighted ALPHA-stable part are + 3.8253 0.2005 + + The reduced controller state dynamics matrix Ac is + 9.1900 0.0000 + 0.0000 -34.5297 + + The reduced controller input/state matrix Bc is + -11.9593 + 86.3137 + + The reduced controller state/output matrix Cc is + 2.8955 -1.3566 + + The reduced controller input/output matrix Dc is + 0.0000 +%} \ No newline at end of file Modified: trunk/octave-forge/extra/control-devel/src/slsb16ad.cc =================================================================== --- trunk/octave-forge/extra/control-devel/src/slsb16ad.cc 2011-11-29 14:00:17 UTC (rev 9220) +++ trunk/octave-forge/extra/control-devel/src/slsb16ad.cc 2011-11-29 14:15:32 UTC (rev 9221) @@ -65,7 +65,7 @@ int nargin = args.length (); octave_value_list retval; - if (nargin != 25) + if (nargin != 19) { print_usage (); } This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2011-11-30 13:05:20
|
Revision: 9224 http://octave.svn.sourceforge.net/octave/?rev=9224&view=rev Author: paramaniac Date: 2011-11-30 13:05:10 +0000 (Wed, 30 Nov 2011) Log Message: ----------- control-devel: quicksave draft code Modified Paths: -------------- trunk/octave-forge/extra/control-devel/inst/bstmodred.m Added Paths: ----------- trunk/octave-forge/extra/control-devel/src/slsb16bd.cc Modified: trunk/octave-forge/extra/control-devel/inst/bstmodred.m =================================================================== --- trunk/octave-forge/extra/control-devel/inst/bstmodred.m 2011-11-30 08:36:23 UTC (rev 9223) +++ trunk/octave-forge/extra/control-devel/inst/bstmodred.m 2011-11-30 13:05:10 UTC (rev 9224) @@ -61,8 +61,8 @@ ## unsolved so far. Nevertheless, balanced truncation and related ## methods can be used to obtain good approximations using this measure. ## -## Available approximation methods are the standard square-root (SR) or -## the accuracy-enhanced balancing-free square-root (BFSR) versions of +## Available approximation methods are the accuracy-enhancing square-root (SR) +## or the balancing-free square-root (BFSR) versions of ## the Balance & Truncate (BTA) or Singular Perturbation Approximation (SPA) ## model reduction methods for the ALPHA-stable part of the system. ## Added: trunk/octave-forge/extra/control-devel/src/slsb16bd.cc =================================================================== --- trunk/octave-forge/extra/control-devel/src/slsb16bd.cc (rev 0) +++ trunk/octave-forge/extra/control-devel/src/slsb16bd.cc 2011-11-30 13:05:10 UTC (rev 9224) @@ -0,0 +1,319 @@ +/* + +Copyright (C) 2011 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/>. + +TODO +Uses SLICOT SB16BD by courtesy of NICONET e.V. +<http://www.slicot.org> + +Author: Lukas Reichlin <luk...@gm...> +Created: November 2011 +Version: 0.1 + +*/ + +#include <octave/oct.h> +#include <f77-fcn.h> +#include "common.cc" + +extern "C" +{ + int F77_FUNC (sb16bd, SB16BD) + (char& DICO, char& JOBD, char& JOBMR, char& JOBCF, + char& EQUIL, char& ORDSEL, + int& N, int& M, int& P, + int& NCR, + double* A, int& LDA, + double* B, int& LDB, + double* C, int& LDC, + double* D, int& LDD, + double* F, int& LDF, + double* G, int& LDG, + double* DC, int& LDDC, + double* HSV, + double& TOL1, double& TOL2, + int* IWORK, + double* DWORK, int& LDWORK, + int& IWARN, int& INFO); +} + +DEFUN_DLD (slsb16bd, args, nargout, + "-*- texinfo -*-\n\ +Slicot SB16BD Release 5.0\n\ +No argument checking.\n\ +For internal use only.") +{ + int nargin = args.length (); + octave_value_list retval; + + if (nargin != 19) + { + print_usage (); + } + else + { + // arguments in + char dico; + char jobc; + char jobo; + char jobmr; + char weight; + char equil; + char ordsel; + + Matrix a = args(0).matrix_value (); + Matrix b = args(1).matrix_value (); + Matrix c = args(2).matrix_value (); + Matrix d = args(3).matrix_value (); + + const int idico = args(4).int_value (); + const int iequil = args(5).int_value (); + int ncr = args(6).int_value (); + const int iordsel = args(7).int_value (); + double alpha = args(8).double_value (); + const int ijobmr = args(9).int_value (); + + Matrix ac = args(10).matrix_value (); + Matrix bc = args(11).matrix_value (); + Matrix cc = args(12).matrix_value (); + Matrix dc = args(13).matrix_value (); + + const int iweight = args(14).int_value (); + const int ijobc = args(15).int_value (); + const int ijobo = args(16).int_value (); + + double tol1 = args(17).double_value (); + double tol2 = args(18).double_value (); + + if (idico == 0) + dico = 'C'; + else + dico = 'D'; + + if (iequil == 0) + equil = 'S'; + else + equil = 'N'; + + if (iordsel == 0) + ordsel = 'F'; + else + ordsel = 'A'; + + if (ijobc == 0) + jobc = 'S'; + else + jobc = 'E'; + + if (ijobo == 0) + jobo = 'S'; + else + jobo = 'E'; + + switch (ijobmr) + { + case 0: + jobmr = 'B'; + break; + case 1: + jobmr = 'F'; + break; + case 2: + jobmr = 'S'; + break; + case 3: + jobmr = 'P'; + break; + default: + error ("slsb16bd: argument jobmr invalid"); + } + + switch (iweight) + { + case 0: + weight = 'N'; + break; + case 1: + weight = 'O'; + break; + case 2: + weight = 'I'; + break; + case 3: + weight = 'P'; + break; + default: + error ("slsb16bd: argument weight invalid"); + } + + int n = a.rows (); // n: number of states + int m = b.columns (); // m: number of inputs + int p = c.rows (); // p: number of outputs + + int nc = ac.rows (); + + int lda = max (1, n); + int ldb = max (1, n); + int ldc = max (1, p); + int ldd = max (1, p); + + int ldac = max (1, nc); + int ldbc = max (1, nc); + int ldcc = max (1, m); + int lddc = max (1, m); + + // arguments out + int ncs; + ColumnVector hsvc (n); + + // workspace + int liwork; + int liwrk1; + int liwrk2; + + switch (jobmr) + { + case 'B': + liwrk1 = 0; + break; + case 'F': + liwrk1 = nc; + break; + default: + liwrk1 = 2*nc; + } + + if (weight == 'N') + liwrk2 = 0; + else + liwrk2 = 2*(m+p); + + liwork = max (1, liwrk1, liwrk2); + + int ldwork; + int lfreq; + int lsqred; + + if (weight == 'N') + { + if (equil == 'N') // if WEIGHT = 'N' and EQUIL = 'N' + lfreq = nc*(max (m, p) + 5); + else // if WEIGHT = 'N' and EQUIL = 'S' + lfreq = max (n, nc*(max (m, p) + 5)); + + } + else // if WEIGHT = 'I' or 'O' or 'P' + { + lfreq = (n+nc)*(n+nc+2*m+2*p) + + max ((n+nc)*(n+nc+max(n+nc,m,p)+7), (m+p)*(m+p+4)); + } + + lsqred = max (1, 2*nc*nc+5*nc); + ldwork = 2*nc*nc + max (1, lfreq, lsqred); + + OCTAVE_LOCAL_BUFFER (int, iwork, liwork); + OCTAVE_LOCAL_BUFFER (double, dwork, ldwork); + + // error indicators + int iwarn = 0; + int info = 0; + + + // SLICOT routine SB16BD + F77_XFCN (sb16bd, SB16BD, + (dico, jobd, jobmr, jobcf, + equil, ordsel, + n, m, p, + ncr, + a.fortran_vec (), lda, + b.fortran_vec (), ldb, + c.fortran_vec (), ldc, + d.fortran_vec (), ldd, + f.fortran_vec (), ldf, + g.fortran_vec (), ldg, + dc.fortran_vec (), lddc, + hsv.fortran_vec (), + tol1, tol2, + iwork, + dwork, ldwork, + iwarn, info)); + + + if (f77_exception_encountered) + error ("conred: exception in SLICOT subroutine SB16BD"); + + if (info != 0) + { + if (info < 0) + error ("conred: the %d-th argument had an invalid value", info); + else + { + switch (info) + { + case 1: + error ("conred: 1: the reduction of A+G*C to a real Schur form " + "failed"); + case 2: + error ("conred: 2: the matrix A+G*C is not stable (if DICO = 'C'), " + "or not convergent (if DICO = 'D')"); + case 3: + error ("conred: 3: the computation of Hankel singular values failed"); + case 4: + error ("conred: 4: the reduction of A+B*F to a real Schur form " + "failed"); + case 5: + error ("conred: 5: the matrix A+B*F is not stable (if DICO = 'C'), " + "or not convergent (if DICO = 'D')"); + default: + error ("conred: unknown error, info = %d", info); + } + } + } + + if (iwarn != 0) + { + switch (iwarn) + { + case 1: + warning ("conred: 1: with ORDSEL = 'F', the selected order NCR is " + "greater than the order of a minimal " + "realization of the controller."); + break; + default: + warning ("conred: unknown warning, iwarn = %d", iwarn); + } + } + + // resize + ac.resize (ncr, ncr); + bc.resize (ncr, p); // p: number of plant outputs + cc.resize (m, ncr); // m: number of plant inputs + hsvc.resize (ncs); + + // return values + retval(0) = ac; + retval(1) = bc; + retval(2) = cc; + retval(3) = dc; + retval(4) = octave_value (ncr); + retval(5) = hsvc; + retval(6) = octave_value (ncs); + } + + return retval; +} This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2011-11-30 14:04:12
|
Revision: 9225 http://octave.svn.sourceforge.net/octave/?rev=9225&view=rev Author: paramaniac Date: 2011-11-30 14:04:01 +0000 (Wed, 30 Nov 2011) Log Message: ----------- control-devel: finish sb16bd Modified Paths: -------------- trunk/octave-forge/extra/control-devel/devel/makefile_conred.m trunk/octave-forge/extra/control-devel/src/slsb16bd.cc Modified: trunk/octave-forge/extra/control-devel/devel/makefile_conred.m =================================================================== --- trunk/octave-forge/extra/control-devel/devel/makefile_conred.m 2011-11-30 13:05:10 UTC (rev 9224) +++ trunk/octave-forge/extra/control-devel/devel/makefile_conred.m 2011-11-30 14:04:01 UTC (rev 9225) @@ -21,7 +21,9 @@ SB03OT.f MB04OY.f MB03QD.f MB04ND.f MB03QY.f \ SB03OR.f SB03OY.f SB04PX.f MB04NY.f SB03OV.f -mkoctfile SB16BD.f AB09AD.f AB09BD.f SB08GD.f SB08HD.f \ +mkoctfile "-Wl,-framework" "-Wl,vecLib" \ + slsb16bd.cc \ + SB16BD.f AB09AD.f AB09BD.f SB08GD.f SB08HD.f \ TB01ID.f AB09AX.f MA02GD.f AB09BX.f TB01WD.f \ MA02DD.f MB03UD.f select.f AB09DD.f SB03OU.f \ MA02AD.f SB03OT.f MB04ND.f MB04OD.f SB03OR.f \ Modified: trunk/octave-forge/extra/control-devel/src/slsb16bd.cc =================================================================== --- trunk/octave-forge/extra/control-devel/src/slsb16bd.cc 2011-11-30 13:05:10 UTC (rev 9224) +++ trunk/octave-forge/extra/control-devel/src/slsb16bd.cc 2011-11-30 14:04:01 UTC (rev 9225) @@ -61,7 +61,7 @@ int nargin = args.length (); octave_value_list retval; - if (nargin != 19) + if (nargin != 15) { print_usage (); } @@ -69,10 +69,9 @@ { // arguments in char dico; - char jobc; - char jobo; + char jobd; char jobmr; - char weight; + char jobcf; char equil; char ordsel; @@ -85,21 +84,17 @@ const int iequil = args(5).int_value (); int ncr = args(6).int_value (); const int iordsel = args(7).int_value (); - double alpha = args(8).double_value (); + const int ijobd = args(8).int_value (); const int ijobmr = args(9).int_value (); - Matrix ac = args(10).matrix_value (); - Matrix bc = args(11).matrix_value (); - Matrix cc = args(12).matrix_value (); - Matrix dc = args(13).matrix_value (); + Matrix f = args(10).matrix_value (); + Matrix g = args(11).matrix_value (); + + const int ijobcf = args(12).int_value (); - const int iweight = args(14).int_value (); - const int ijobc = args(15).int_value (); - const int ijobo = args(16).int_value (); + double tol1 = args(13).double_value (); + double tol2 = args(14).double_value (); - double tol1 = args(17).double_value (); - double tol2 = args(18).double_value (); - if (idico == 0) dico = 'C'; else @@ -115,15 +110,15 @@ else ordsel = 'A'; - if (ijobc == 0) - jobc = 'S'; + if (ijobd == 0) + jobd = 'Z'; else - jobc = 'E'; + jobd = 'D'; - if (ijobo == 0) - jobo = 'S'; + if (ijobcf == 0) + jobcf = 'L'; else - jobo = 'E'; + jobcf = 'R'; switch (ijobmr) { @@ -143,88 +138,64 @@ error ("slsb16bd: argument jobmr invalid"); } - switch (iweight) - { - case 0: - weight = 'N'; - break; - case 1: - weight = 'O'; - break; - case 2: - weight = 'I'; - break; - case 3: - weight = 'P'; - break; - default: - error ("slsb16bd: argument weight invalid"); - } int n = a.rows (); // n: number of states int m = b.columns (); // m: number of inputs int p = c.rows (); // p: number of outputs - - int nc = ac.rows (); int lda = max (1, n); int ldb = max (1, n); int ldc = max (1, p); - int ldd = max (1, p); + int ldd; + + if (jobd == 'Z') + ldd = 1; + else + ldd = max (1, p); - int ldac = max (1, nc); - int ldbc = max (1, nc); - int ldcc = max (1, m); + int ldf = max (1, m); + int ldg = max (1, n); int lddc = max (1, m); // arguments out - int ncs; - ColumnVector hsvc (n); + Matrix dc (lddc, p); + ColumnVector hsv (n); // workspace int liwork; - int liwrk1; - int liwrk2; + int pm; + int ldwork; + int lwr = max (1, n*(2*n+max(n,m+p)+5)+n*(n+1)/2); + switch (jobmr) { case 'B': - liwrk1 = 0; + pm = 0; break; case 'F': - liwrk1 = nc; + pm = n; break; - default: - liwrk1 = 2*nc; + default: // if JOBMR = 'S' or 'P' + pm = max (1, 2*n); } - if (weight == 'N') - liwrk2 = 0; - else - liwrk2 = 2*(m+p); - - liwork = max (1, liwrk1, liwrk2); - - int ldwork; - int lfreq; - int lsqred; - - if (weight == 'N') + if (ordsel == 'F' && ncr == n) { - if (equil == 'N') // if WEIGHT = 'N' and EQUIL = 'N' - lfreq = nc*(max (m, p) + 5); - else // if WEIGHT = 'N' and EQUIL = 'S' - lfreq = max (n, nc*(max (m, p) + 5)); - + liwork = 0; + ldwork = p*n; } - else // if WEIGHT = 'I' or 'O' or 'P' + else if (jobcf == 'L') { - lfreq = (n+nc)*(n+nc+2*m+2*p) + - max ((n+nc)*(n+nc+max(n+nc,m,p)+7), (m+p)*(m+p+4)); + liwork = max (pm, m); + ldwork = (n+m)*(m+p) + max (lwr, 4*m); } + else // if JOBCF = 'R' + { + liwork = max (pm, p); + ldwork = (n+p)*(m+p) + max (lwr, 4*p); + } - lsqred = max (1, 2*nc*nc+5*nc); - ldwork = 2*nc*nc + max (1, lfreq, lsqred); OCTAVE_LOCAL_BUFFER (int, iwork, liwork); OCTAVE_LOCAL_BUFFER (double, dwork, ldwork); @@ -300,19 +271,18 @@ } // resize - ac.resize (ncr, ncr); - bc.resize (ncr, p); // p: number of plant outputs - cc.resize (m, ncr); // m: number of plant inputs - hsvc.resize (ncs); + a.resize (ncr, ncr); // Ac + g.resize (ncr, p); // Bc + f.resize (m, ncr); // Cc + dc.resize (m, p); // Dc // return values - retval(0) = ac; - retval(1) = bc; - retval(2) = cc; + retval(0) = a; + retval(1) = g; + retval(2) = f; retval(3) = dc; retval(4) = octave_value (ncr); - retval(5) = hsvc; - retval(6) = octave_value (ncs); + retval(5) = hsv; } return retval; This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2011-11-30 14:54:26
|
Revision: 9227 http://octave.svn.sourceforge.net/octave/?rev=9227&view=rev Author: paramaniac Date: 2011-11-30 14:54:12 +0000 (Wed, 30 Nov 2011) Log Message: ----------- control-devel: add slicot sb16cd Modified Paths: -------------- trunk/octave-forge/extra/control-devel/devel/makefile_conred.m Added Paths: ----------- trunk/octave-forge/extra/control-devel/src/slsb16cd.cc Modified: trunk/octave-forge/extra/control-devel/devel/makefile_conred.m =================================================================== --- trunk/octave-forge/extra/control-devel/devel/makefile_conred.m 2011-11-30 14:18:03 UTC (rev 9226) +++ trunk/octave-forge/extra/control-devel/devel/makefile_conred.m 2011-11-30 14:54:12 UTC (rev 9227) @@ -29,7 +29,9 @@ MA02AD.f SB03OT.f MB04ND.f MB04OD.f SB03OR.f \ SB03OY.f SB04PX.f MB04NY.f MB04OY.f SB03OV.f -mkoctfile SB16CD.f SB16CY.f AB09IX.f SB03OD.f MB02UD.f \ +mkoctfile "-Wl,-framework" "-Wl,vecLib" \ + slsb16cd.cc \ + SB16CD.f SB16CY.f AB09IX.f SB03OD.f MB02UD.f \ AB09DD.f MA02AD.f MB03UD.f select.f SB03OU.f \ MB01SD.f SB03OT.f MB04ND.f MB04OD.f SB03OR.f \ SB03OY.f SB04PX.f MB04NY.f MB04OY.f SB03OV.f Added: trunk/octave-forge/extra/control-devel/src/slsb16cd.cc =================================================================== --- trunk/octave-forge/extra/control-devel/src/slsb16cd.cc (rev 0) +++ trunk/octave-forge/extra/control-devel/src/slsb16cd.cc 2011-11-30 14:54:12 UTC (rev 9227) @@ -0,0 +1,261 @@ +/* + +Copyright (C) 2011 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/>. + +TODO +Uses SLICOT SB16CD by courtesy of NICONET e.V. +<http://www.slicot.org> + +Author: Lukas Reichlin <luk...@gm...> +Created: November 2011 +Version: 0.1 + +*/ + +#include <octave/oct.h> +#include <f77-fcn.h> +#include "common.cc" + +extern "C" +{ + int F77_FUNC (sb16cd, SB16CD) + (char& DICO, char& JOBD, char& JOBMR, char& JOBCF, + char& ORDSEL, + int& N, int& M, int& P, + int& NCR, + double* A, int& LDA, + double* B, int& LDB, + double* C, int& LDC, + double* D, int& LDD, + double* F, int& LDF, + double* G, int& LDG, + double* HSV, + double& TOL, + int* IWORK, + double* DWORK, int& LDWORK, + int& IWARN, int& INFO); +} + +DEFUN_DLD (slsb16cd, args, nargout, + "-*- texinfo -*-\n\ +Slicot SB16CD Release 5.0\n\ +No argument checking.\n\ +For internal use only.") +{ + int nargin = args.length (); + octave_value_list retval; + + if (nargin != 15) + { + print_usage (); + } + else + { + // arguments in + char dico; + char jobd; + char jobmr; + char jobcf; + char ordsel; + + Matrix a = args(0).matrix_value (); + Matrix b = args(1).matrix_value (); + Matrix c = args(2).matrix_value (); + Matrix d = args(3).matrix_value (); + + const int idico = args(4).int_value (); + int ncr = args(5).int_value (); + const int iordsel = args(6).int_value (); + const int ijobd = args(7).int_value (); + const int ijobmr = args(8).int_value (); + + Matrix f = args(9).matrix_value (); + Matrix g = args(10).matrix_value (); + + const int ijobcf = args(11).int_value (); + double tol = args(12).double_value (); + + if (idico == 0) + dico = 'C'; + else + dico = 'D'; + + if (iordsel == 0) + ordsel = 'F'; + else + ordsel = 'A'; + + if (ijobd == 0) + jobd = 'Z'; + else + jobd = 'D'; + + if (ijobcf == 0) + jobcf = 'L'; + else + jobcf = 'R'; + + switch (ijobmr) + { + case 0: + jobmr = 'B'; + break; + case 1: + jobmr = 'F'; + break; + default: + error ("slsb16cd: argument jobmr invalid"); + } + + + int n = a.rows (); // n: number of states + int m = b.columns (); // m: number of inputs + int p = c.rows (); // p: number of outputs + + int lda = max (1, n); + int ldb = max (1, n); + int ldc = max (1, p); + int ldd; + + if (jobd == 'Z') + ldd = 1; + else + ldd = max (1, p); + + int ldf = max (1, m); + int ldg = max (1, n); + + // arguments out + ColumnVector hsv (n); + + // workspace + int liwork; + + if (jobmr == 'B') + liwork = 0; + else // if JOBMR = 'F' + liwork = n; + + int ldwork; + int mp; + + if (jobcf == 'L') + mp = m; + else // if JOBCF = 'R' + mp = p; + + ldwork = 2*n*n + max (1, 2*n*n + 5*n, n*max(m,p), + n*(n + max(n,mp) + min(n,mp) + 6)); + + + OCTAVE_LOCAL_BUFFER (int, iwork, liwork); + OCTAVE_LOCAL_BUFFER (double, dwork, ldwork); + + // error indicators + int iwarn = 0; + int info = 0; + + + // SLICOT routine SB16CD + F77_XFCN (sb16cd, SB16CD, + (dico, jobd, jobmr, jobcf, + ordsel, + n, m, p, + ncr, + a.fortran_vec (), lda, + b.fortran_vec (), ldb, + c.fortran_vec (), ldc, + d.fortran_vec (), ldd, + f.fortran_vec (), ldf, + g.fortran_vec (), ldg, + hsv.fortran_vec (), + tol, + iwork, + dwork, ldwork, + iwarn, info)); + + + if (f77_exception_encountered) + error ("conred: exception in SLICOT subroutine SB16CD"); + + if (info != 0) + { + if (info < 0) + error ("conred: the %d-th argument had an invalid value", info); + else + { + switch (info) + { + case 1: + error ("conred: 1: eigenvalue computation failure"); + case 2: + error ("conred: 2: the matrix A+G*C is not stable"); + case 3: + error ("conred: 3: the matrix A+B*F is not stable"); + case 4: + error ("conred: 4: the Lyapunov equation for computing the " + "observability Grammian is (nearly) singular"); + case 5: + error ("conred: 5: the Lyapunov equation for computing the " + "controllability Grammian is (nearly) singular"); + case 6: + error ("conred: 6: the computation of Hankel singular values failed"); + default: + error ("conred: unknown error, info = %d", info); + } + } + } + + if (iwarn != 0) + { + switch (iwarn) + { + case 1: + warning ("conred: 1: with ORDSEL = 'F', the selected order NCR is " + "greater than the order of a minimal realization " + "of the controller."); + break; + case 2: + warning ("conred: 2: with ORDSEL = 'F', the selected order NCR " + "corresponds to repeated singular values, which are " + "neither all included nor all excluded from the " + "reduced controller. In this case, the resulting NCR " + "is set automatically to the largest value such that " + "HSV(NCR) > HSV(NCR+1)."); + break; + default: + warning ("conred: unknown warning, iwarn = %d", iwarn); + } + } + + // resize + a.resize (ncr, ncr); // Ac + g.resize (ncr, p); // Bc + f.resize (m, ncr); // Cc + // Dc = 0 + + // return values + retval(0) = a; + retval(1) = g; + retval(2) = f; + retval(3) = octave_value (ncr); + retval(4) = hsv; + } + + return retval; +} This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2011-11-30 15:04:00
|
Revision: 9228 http://octave.svn.sourceforge.net/octave/?rev=9228&view=rev Author: paramaniac Date: 2011-11-30 15:03:54 +0000 (Wed, 30 Nov 2011) Log Message: ----------- control-devel: get slicot sb16cd working Modified Paths: -------------- trunk/octave-forge/extra/control-devel/src/slsb16cd.cc Added Paths: ----------- trunk/octave-forge/extra/control-devel/devel/test_sb16cd.m Added: trunk/octave-forge/extra/control-devel/devel/test_sb16cd.m =================================================================== --- trunk/octave-forge/extra/control-devel/devel/test_sb16cd.m (rev 0) +++ trunk/octave-forge/extra/control-devel/devel/test_sb16cd.m 2011-11-30 15:03:54 UTC (rev 9228) @@ -0,0 +1,108 @@ +%{ + SB16CD EXAMPLE PROGRAM DATA (Continuous system) + 8 1 1 2 0.1E0 C D F R F + + N, M, P, NCR, TOL, + DICO, JOBD, JOBMR, JOBCF, ORDSEL + + 0 1.0000 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 + 0 0 -0.0150 0.7650 0 0 0 0 + 0 0 -0.7650 -0.0150 0 0 0 0 + 0 0 0 0 -0.0280 1.4100 0 0 + 0 0 0 0 -1.4100 -0.0280 0 0 + 0 0 0 0 0 0 -0.0400 1.850 + 0 0 0 0 0 0 -1.8500 -0.040 + 0.0260 + -0.2510 + 0.0330 + -0.8860 + -4.0170 + 0.1450 + 3.6040 + 0.2800 + -.996 -.105 0.261 .009 -.001 -.043 0.002 -0.026 + 0.0 +4.472135954999638e-002 6.610515358414598e-001 4.698598960657579e-003 3.601363251422058e-001 1.032530880771415e-001 -3.754055214487997e-002 -4.268536964759344e-002 3.287284547842979e-002 + 4.108939884667451e-001 + 8.684600000000012e-002 + 3.852317308197148e-004 + -3.619366874815911e-003 + -8.803722876359955e-003 + 8.420521094001852e-003 + 1.234944428038507e-003 + 4.263205617645322e-003 +%} + +a = [ 0 1.0000 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 + 0 0 -0.0150 0.7650 0 0 0 0 + 0 0 -0.7650 -0.0150 0 0 0 0 + 0 0 0 0 -0.0280 1.4100 0 0 + 0 0 0 0 -1.4100 -0.0280 0 0 + 0 0 0 0 0 0 -0.0400 1.850 + 0 0 0 0 0 0 -1.8500 -0.040 ]; + +b = [ 0.0260 + -0.2510 + 0.0330 + -0.8860 + -4.0170 + 0.1450 + 3.6040 + 0.2800 ]; + +c = [ -.996 -.105 0.261 .009 -.001 -.043 0.002 -0.026 ]; + +d = [ 0.0 ]; + +f = [ 4.472135954999638e-002 6.610515358414598e-001 4.698598960657579e-003 3.601363251422058e-001 1.032530880771415e-001 -3.754055214487997e-002 -4.268536964759344e-002 3.287284547842979e-002 ]; + +g = [ 4.108939884667451e-001 + 8.684600000000012e-002 + 3.852317308197148e-004 + -3.619366874815911e-003 + -8.803722876359955e-003 + 8.420521094001852e-003 + 1.234944428038507e-003 + 4.263205617645322e-003 ]; + +%{ + 8 1 1 2 0.1E0 C D F R F + + N, M, P, NCR, TOL, + DICO, JOBD, JOBMR, JOBCF, ORDSEL +%} + + +tol = 0.1; # tol1 +dico = 0; +jobd = 1; +jobmr = 1; +jobcf = 1; +ordsel = 0; +ncr = 2; + + +[ac, bc, cc] = slsb16cd (a, b, c, d, dico, ncr, ordsel, jobd, jobmr, \ + f, g, jobcf, tol) + +%{ + SB16CD EXAMPLE PROGRAM RESULTS + + The order of reduced controller = 2 + + The frequency-weighted Hankel singular values are: + 3.3073 0.7274 0.1124 0.0784 0.0242 0.0182 0.0101 0.0094 + + The reduced controller state dynamics matrix Ac is + -0.4334 0.4884 + -0.1950 -0.1093 + + The reduced controller input/state matrix Bc is + -0.4231 + -0.1785 + + The reduced controller state/output matrix Cc is + -0.0326 -0.2307 +%} \ No newline at end of file Modified: trunk/octave-forge/extra/control-devel/src/slsb16cd.cc =================================================================== --- trunk/octave-forge/extra/control-devel/src/slsb16cd.cc 2011-11-30 14:54:12 UTC (rev 9227) +++ trunk/octave-forge/extra/control-devel/src/slsb16cd.cc 2011-11-30 15:03:54 UTC (rev 9228) @@ -60,7 +60,7 @@ int nargin = args.length (); octave_value_list retval; - if (nargin != 15) + if (nargin != 13) { print_usage (); } This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2011-11-30 19:37:51
|
Revision: 9229 http://octave.svn.sourceforge.net/octave/?rev=9229&view=rev Author: paramaniac Date: 2011-11-30 19:37:41 +0000 (Wed, 30 Nov 2011) Log Message: ----------- control-devel: improve argument handling for bstmodred Modified Paths: -------------- trunk/octave-forge/extra/control-devel/devel/test_bstmodred.m trunk/octave-forge/extra/control-devel/inst/__opt2cell__.m trunk/octave-forge/extra/control-devel/inst/bstmodred.m Modified: trunk/octave-forge/extra/control-devel/devel/test_bstmodred.m =================================================================== --- trunk/octave-forge/extra/control-devel/devel/test_bstmodred.m 2011-11-30 15:03:54 UTC (rev 9228) +++ trunk/octave-forge/extra/control-devel/devel/test_bstmodred.m 2011-11-30 19:37:41 UTC (rev 9229) @@ -26,3 +26,12 @@ sysr = bstmodred (sys, "beta", 1.0, "tol1", 0.1, "tol2", 0.0) [Ao, Bo, Co, Do] = ssdata (sysr); + + +opt = options ("beta", 1.0, "tol1", 0.1, "tol2", 0.0) +sysr = bstmodred (sys, opt) + + +sysr = bstmodred (sys, 5, "beta", 1.0, "tol1", 0.1, "tol2", 0.0) + +sysr = bstmodred (sys, 5, opt) Modified: trunk/octave-forge/extra/control-devel/inst/__opt2cell__.m =================================================================== --- trunk/octave-forge/extra/control-devel/inst/__opt2cell__.m 2011-11-30 15:03:54 UTC (rev 9228) +++ trunk/octave-forge/extra/control-devel/inst/__opt2cell__.m 2011-11-30 19:37:41 UTC (rev 9229) @@ -32,6 +32,6 @@ key = fieldnames (opt); val = struct2cell (opt); - c = [key.'; val.'](:); # reshape to {key1; val1; key2; val2; ...} + c = [key.'; val.'](:).'; # reshape to {key1, val1, key2, val2, ...} endfunction \ No newline at end of file Modified: trunk/octave-forge/extra/control-devel/inst/bstmodred.m =================================================================== --- trunk/octave-forge/extra/control-devel/inst/bstmodred.m 2011-11-30 15:03:54 UTC (rev 9228) +++ trunk/octave-forge/extra/control-devel/inst/bstmodred.m 2011-11-30 19:37:41 UTC (rev 9229) @@ -16,9 +16,11 @@ ## along with LTI Syncope. If not, see <http://www.gnu.org/licenses/>. ## -*- texinfo -*- -## @deftypefn{Function File} {[@var{Gr}, @var{info}] =} bstmodred (@var{G}) -## @deftypefnx{Function File} {[@var{Gr}, @var{info}] =} bstmodred (@var{G}, @dots{}) -## @deftypefnx{Function File} {[@var{Gr}, @var{info}] =} bstmodred (@var{G}, @var{opt}) +## @deftypefn{Function File} {[@var{Gr}, @var{info}] =} bstmodred (@var{G}, @dots{}) +## @deftypefnx{Function File} {[@var{Gr}, @var{info}] =} bstmodred (@var{G}, @var{nr}, @dots{}) +## @deftypefnx{Function File} {[@var{Gr}, @var{info}] =} bstmodred (@var{G}, @var{opt}, @dots{}) +## @deftypefnx{Function File} {[@var{Gr}, @var{info}] =} bstmodred (@var{G}, @var{nr}, @var{opt}, @dots{}) +## ## Model order reduction by Balanced Stochastic Truncation (BST) method. ## The aim of model reduction is to find an LTI system @var{Gr} of order ## @var{nr} (nr << n) such that the input-output behaviour of @var{Gr} @@ -109,9 +111,9 @@ ## @item G ## LTI model to be reduced. ## @item @dots{} -## Pairs of keys and values. +## Optional pairs of keys and values. ## @item opt -## Struct with keys as field names. +## Optional struct with keys as field names. ## @end table ## ## @strong{Outputs} @@ -228,11 +230,24 @@ error ("bstmodred: first argument must be an LTI system"); endif - if (nargin == 2) - varargin = __opt2cell__ (varargin{1}); - elseif (rem (nargin-1, 2)) - error ("bstmodred: properties and values must come in pairs"); + if (nargin > 1) # bstmodred (G, ...) + if (isstruct (varargin{1})) # bstmodred (G, opt, ...) + varargin = horzcat (__opt2cell__ (varargin{1}), varargin(2:end)); + elseif (is_real_scalar (varargin{1})) # bstmodred (G, nr) + varargin = horzcat ({"order"}, varargin); + if (nargin > 2 && isstruct (varargin{3})) # bstmodred (G, nr, ...) + varargin = horzcat (__opt2cell__ (varargin{3}), varargin([1:2, 4:end])); + ## varargin(1:2) placed after opt such that nr from bstmodred (G, nr, opt, ...) + ## overrides possible nr's inside opt struct (later keys override former keys) + endif + endif endif + + narg = numel (varargin); + + if (rem (narg, 2)) + error ("bstmodred: keys and values must come in pairs"); + endif [a, b, c, d, tsam, scaled] = ssdata (sys); dt = isdt (sys); @@ -247,11 +262,11 @@ job = 1; ## handle properties and values - for k = 1 : 2 : (nargin-1) + for k = 1 : 2 : narg prop = lower (varargin{k}); val = varargin{k+1}; switch (prop) - case {"order", "n", "nr"} + case {"order", "nr"} [nr, ordsel] = __modred_check_order__ (val); case "tol1" @@ -360,4 +375,3 @@ %! Me = [Ae, Be; Ce, De]; %! %!assert (Mo, Me, 1e-4); - This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |