1 /* 2 3 Copyright (C) 2009-2016 Lukas F. Reichlin 4 5 This file is part of LTI Syncope. 6 7 LTI Syncope is free software: you can redistribute it and/or modify 8 it under the terms of the GNU General Public License as published by 9 the Free Software Foundation, either version 3 of the License, or 10 (at your option) any later version. 11 12 LTI Syncope is distributed in the hope that it will be useful, 13 but WITHOUT ANY WARRANTY; without even the implied warranty of 14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 15 GNU General Public License for more details. 16 17 You should have received a copy of the GNU General Public License 18 along with LTI Syncope. If not, see <http://www.gnu.org/licenses/>. 19 20 Controller reduction based on Balance & Truncate (B&T) or 21 Singular Perturbation Approximation (SPA) method. 22 Uses SLICOT SB16AD by courtesy of NICONET e.V. 23 <http://www.slicot.org> 24 25 Author: Lukas Reichlin <lukas.reichlin@gmail.com> 26 Created: November 2011 27 Version: 0.2 28 29 */ 30 31 #include <octave/oct.h> 32 #include "common.h" 33 34 extern "C" 35 { 36 int F77_FUNC (sb16ad, SB16AD) 37 (char& DICO, char& JOBC, char& JOBO, char& JOBMR, 38 char& WEIGHT, char& EQUIL, char& ORDSEL, 39 F77_INT& N, F77_INT& M, F77_INT& P, 40 F77_INT& NC, F77_INT& NCR, 41 double& ALPHA, 42 double* A, F77_INT& LDA, 43 double* B, F77_INT& LDB, 44 double* C, F77_INT& LDC, 45 double* D, F77_INT& LDD, 46 double* AC, F77_INT& LDAC, 47 double* BC, F77_INT& LDBC, 48 double* CC, F77_INT& LDCC, 49 double* DC, F77_INT& LDDC, 50 F77_INT& NCS, 51 double* HSVC, 52 double& TOL1, double& TOL2, 53 F77_INT* IWORK, 54 double* DWORK, F77_INT& LDWORK, 55 F77_INT& IWARN, F77_INT& INFO); 56 } 57 58 // PKG_ADD: autoload ("__sl_sb16ad__", "__control_slicot_functions__.oct"); 59 DEFUN_DLD (__sl_sb16ad__, args, nargout, 60 "-*- texinfo -*-\n\ 61 Slicot SB16AD Release 5.0\n\ 62 No argument checking.\n\ 63 For internal use only.") 64 { 65 octave_idx_type nargin = args.length (); 66 octave_value_list retval; 67 68 if (nargin != 19) 69 { 70 print_usage (); 71 } 72 else 73 { 74 // arguments in 75 char dico; 76 char jobc; 77 char jobo; 78 char jobmr; 79 char weight; 80 char equil; 81 char ordsel; 82 83 Matrix a = args(0).matrix_value (); 84 Matrix b = args(1).matrix_value (); 85 Matrix c = args(2).matrix_value (); 86 Matrix d = args(3).matrix_value (); 87 88 const F77_INT idico = args(4).int_value (); 89 const F77_INT iequil = args(5).int_value (); 90 F77_INT ncr = args(6).int_value (); 91 const F77_INT iordsel = args(7).int_value (); 92 double alpha = args(8).double_value (); 93 const F77_INT ijobmr = args(9).int_value (); 94 95 Matrix ac = args(10).matrix_value (); 96 Matrix bc = args(11).matrix_value (); 97 Matrix cc = args(12).matrix_value (); 98 Matrix dc = args(13).matrix_value (); 99 100 const F77_INT iweight = args(14).int_value (); 101 const F77_INT ijobc = args(15).int_value (); 102 const F77_INT ijobo = args(16).int_value (); 103 104 double tol1 = args(17).double_value (); 105 double tol2 = args(18).double_value (); 106 107 if (idico == 0) 108 dico = 'C'; 109 else 110 dico = 'D'; 111 112 if (iequil == 0) 113 equil = 'S'; 114 else 115 equil = 'N'; 116 117 if (iordsel == 0) 118 ordsel = 'F'; 119 else 120 ordsel = 'A'; 121 122 if (ijobc == 0) 123 jobc = 'S'; 124 else 125 jobc = 'E'; 126 127 if (ijobo == 0) 128 jobo = 'S'; 129 else 130 jobo = 'E'; 131 132 switch (ijobmr) 133 { 134 case 0: 135 jobmr = 'B'; 136 break; 137 case 1: 138 jobmr = 'F'; 139 break; 140 case 2: 141 jobmr = 'S'; 142 break; 143 case 3: 144 jobmr = 'P'; 145 break; 146 default: 147 error ("__sl_sb16ad__: argument jobmr invalid"); 148 } 149 150 switch (iweight) 151 { 152 case 0: 153 weight = 'N'; 154 break; 155 case 1: 156 weight = 'O'; 157 break; 158 case 2: 159 weight = 'I'; 160 break; 161 case 3: 162 weight = 'P'; 163 break; 164 default: 165 error ("__sl_sb16ad__: argument weight invalid"); 166 } 167 168 F77_INT n = TO_F77_INT (a.rows ()); // n: number of states 169 F77_INT m = TO_F77_INT (b.columns ()); // m: number of inputs 170 F77_INT p = TO_F77_INT (c.rows ()); // p: number of outputs 171 172 F77_INT nc = TO_F77_INT (ac.rows ()); 173 174 F77_INT lda = max (1, n); 175 F77_INT ldb = max (1, n); 176 F77_INT ldc = max (1, p); 177 F77_INT ldd = max (1, p); 178 179 F77_INT ldac = max (1, nc); 180 F77_INT ldbc = max (1, nc); 181 F77_INT ldcc = max (1, m); 182 F77_INT lddc = max (1, m); 183 184 // arguments out 185 F77_INT ncs; 186 ColumnVector hsvc (n); 187 188 // workspace 189 F77_INT liwork; 190 F77_INT liwrk1; 191 F77_INT liwrk2; 192 193 switch (jobmr) 194 { 195 case 'B': 196 liwrk1 = 0; 197 break; 198 case 'F': 199 liwrk1 = nc; 200 break; 201 default: 202 liwrk1 = 2*nc; 203 } 204 205 if (weight == 'N') 206 liwrk2 = 0; 207 else 208 liwrk2 = 2*(m+p); 209 210 liwork = max (1, liwrk1, liwrk2); 211 212 F77_INT ldwork; 213 F77_INT lfreq; 214 F77_INT lsqred; 215 216 if (weight == 'N') 217 { 218 if (equil == 'N') // if WEIGHT = 'N' and EQUIL = 'N' 219 lfreq = nc*(max (m, p) + 5); 220 else // if WEIGHT = 'N' and EQUIL = 'S' 221 lfreq = max (n, nc*(max (m, p) + 5)); 222 223 } 224 else // if WEIGHT = 'I' or 'O' or 'P' 225 { 226 lfreq = (n+nc)*(n+nc+2*m+2*p) + 227 max ((n+nc)*(n+nc+max(n+nc,m,p)+7), (m+p)*(m+p+4)); 228 } 229 230 lsqred = max (1, 2*nc*nc+5*nc); 231 ldwork = 2*nc*nc + max (1, lfreq, lsqred); 232 233 OCTAVE_LOCAL_BUFFER (F77_INT, iwork, liwork); 234 OCTAVE_LOCAL_BUFFER (double, dwork, ldwork); 235 236 // error indicators 237 F77_INT iwarn = 0; 238 F77_INT info = 0; 239 240 241 // SLICOT routine SB16AD 242 F77_XFCN (sb16ad, SB16AD, 243 (dico, jobc, jobo, jobmr, 244 weight, equil, ordsel, 245 n, m, p, 246 nc, ncr, 247 alpha, 248 a.fortran_vec (), lda, 249 b.fortran_vec (), ldb, 250 c.fortran_vec (), ldc, 251 d.fortran_vec (), ldd, 252 ac.fortran_vec (), ldac, 253 bc.fortran_vec (), ldbc, 254 cc.fortran_vec (), ldcc, 255 dc.fortran_vec (), lddc, 256 ncs, 257 hsvc.fortran_vec (), 258 tol1, tol2, 259 iwork, 260 dwork, ldwork, 261 iwarn, info)); 262 263 264 if (f77_exception_encountered) 265 error ("conred: exception in SLICOT subroutine SB16AD"); 266 267 static const char* err_msg[] = { 268 "0: OK", 269 "1: the closed-loop system is not well-posed; " 270 "its feedthrough matrix is (numerically) singular", 271 "2: the computation of the real Schur form of the " 272 "closed-loop state matrix failed", 273 "3: the closed-loop state matrix is not stable", 274 "4: the solution of a symmetric eigenproblem failed", 275 "5: the computation of the ordered real Schur form " 276 "of Ac failed", 277 "6: the separation of the ALPHA-stable/unstable " 278 "diagonal blocks failed because of very close eigenvalues", 279 "7: the computation of Hankel singular values failed"}; 280 281 static const char* warn_msg[] = { 282 "0: OK", 283 "1: with ORDSEL = 'F', the selected order NCR is greater " 284 "than NSMIN, the sum of the order of the " 285 "ALPHA-unstable part and the order of a minimal " 286 "realization of the ALPHA-stable part of the given " 287 "controller; in this case, the resulting NCR is set " 288 "equal to NSMIN.", 289 "2: with ORDSEL = 'F', the selected order NCR " 290 "corresponds to repeated singular values for the " 291 "ALPHA-stable part of the controller, which are " 292 "neither all included nor all excluded from the " 293 "reduced model; in this case, the resulting NCR is " 294 "automatically decreased to exclude all repeated " 295 "singular values.", 296 "3: with ORDSEL = 'F', the selected order NCR is less " 297 "than the order of the ALPHA-unstable part of the " 298 "given controller. In this case NCR is set equal to " 299 "the order of the ALPHA-unstable part."}; 300 301 error_msg ("conred", info, 7, err_msg); 302 warning_msg ("conred", iwarn, 3, warn_msg); 303 304 305 // resize 306 ac.resize (ncr, ncr); 307 bc.resize (ncr, p); // p: number of plant outputs 308 cc.resize (m, ncr); // m: number of plant inputs 309 hsvc.resize (ncs); 310 311 // return values 312 retval(0) = ac; 313 retval(1) = bc; 314 retval(2) = cc; 315 retval(3) = dc; 316 retval(4) = octave_value (ncr); 317 retval(5) = hsvc; 318 retval(6) = octave_value (ncs); 319 } 320 321 return retval; 322 } 323