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 Positive feedback controller for a discrete-time system (D == 0). 21 Uses SLICOT SB10KD by courtesy of NICONET e.V. 22 <http://www.slicot.org> 23 24 Author: Lukas Reichlin <lukas.reichlin@gmail.com> 25 Created: July 2011 26 Version: 0.4 27 28 */ 29 30 #include <octave/oct.h> 31 #include "common.h" 32 33 extern "C" 34 { 35 int F77_FUNC (sb10kd, SB10KD) 36 (F77_INT& N, F77_INT& M, F77_INT& NP, 37 double* A, F77_INT& LDA, 38 double* B, F77_INT& LDB, 39 double* C, F77_INT& LDC, 40 double& FACTOR, 41 double* AK, F77_INT& LDAK, 42 double* BK, F77_INT& LDBK, 43 double* CK, F77_INT& LDCK, 44 double* DK, F77_INT& LDDK, 45 double* RCOND, 46 F77_INT* IWORK, 47 double* DWORK, F77_INT& LDWORK, 48 F77_LOGICAL* BWORK, 49 F77_INT& INFO); 50 } 51 52 // PKG_ADD: autoload ("__sl_sb10kd__", "__control_slicot_functions__.oct"); 53 DEFUN_DLD (__sl_sb10kd__, args, nargout, 54 "-*- texinfo -*-\n\ 55 Slicot SB10KD Release 5.0\n\ 56 No argument checking.\n\ 57 For internal use only.") 58 { 59 octave_idx_type nargin = args.length (); 60 octave_value_list retval; 61 62 if (nargin != 4) 63 { 64 print_usage (); 65 } 66 else 67 { 68 // arguments in 69 Matrix a = args(0).matrix_value (); 70 Matrix b = args(1).matrix_value (); 71 Matrix c = args(2).matrix_value (); 72 73 double factor = args(3).double_value (); 74 75 F77_INT n = TO_F77_INT (a.rows ()); // n: number of states 76 F77_INT m = TO_F77_INT (b.columns ()); // m: number of inputs 77 F77_INT np = TO_F77_INT (c.rows ()); // np: number of outputs 78 79 F77_INT lda = max (1, n); 80 F77_INT ldb = max (1, n); 81 F77_INT ldc = max (1, np); 82 83 F77_INT ldak = max (1, n); 84 F77_INT ldbk = max (1, n); 85 F77_INT ldck = max (1, m); 86 F77_INT lddk = max (1, m); 87 88 // arguments out 89 Matrix ak (ldak, n); 90 Matrix bk (ldbk, np); 91 Matrix ck (ldck, n); 92 Matrix dk (lddk, np); 93 ColumnVector rcond (4); 94 95 // workspace 96 F77_INT liwork = 2 * max (n, np+m); 97 F77_INT ldwork = 15*n*n + 6*n + 98 max (14*n+23, 16*n, 2*n+np+m, 3*(np+m)) + 99 max (n*n, 11*n*np + 2*m*m + 8*np*np + 8*m*n + 4*m*np + np); 100 101 OCTAVE_LOCAL_BUFFER (F77_INT, iwork, liwork); 102 OCTAVE_LOCAL_BUFFER (double, dwork, ldwork); 103 OCTAVE_LOCAL_BUFFER (F77_LOGICAL, bwork, 2*n); 104 105 // error indicator 106 F77_INT info; 107 108 109 // SLICOT routine SB10KD 110 F77_XFCN (sb10kd, SB10KD, 111 (n, m, np, 112 a.fortran_vec (), lda, 113 b.fortran_vec (), ldb, 114 c.fortran_vec (), ldc, 115 factor, 116 ak.fortran_vec (), ldak, 117 bk.fortran_vec (), ldbk, 118 ck.fortran_vec (), ldck, 119 dk.fortran_vec (), lddk, 120 rcond.fortran_vec (), 121 iwork, 122 dwork, ldwork, 123 bwork, 124 info)); 125 126 if (f77_exception_encountered) 127 error ("ncfsyn: slsb10kd: exception in SLICOT subroutine SB10KD"); 128 129 static const char* err_msg[] = { 130 "0: OK", 131 "1: the P-Riccati equation is not solved successfully", 132 "2: the Q-Riccati equation is not solved successfully", 133 "3: the X-Riccati equation is not solved successfully", 134 "4: the iteration to compute eigenvalues failed to " 135 "converge", 136 "5: the matrix Rx + Bx'*X*Bx is singular", 137 "6: the closed-loop system is unstable"}; 138 139 error_msg ("ncfsyn", info, 6, err_msg); 140 141 142 // resizing 143 ak.resize (n, n); 144 bk.resize (n, np); 145 ck.resize (m, n); 146 dk.resize (m, np); 147 148 // return values 149 retval(0) = ak; 150 retval(1) = bk; 151 retval(2) = ck; 152 retval(3) = dk; 153 retval(4) = rcond; 154 } 155 156 return retval; 157 } 158