octave_packages/control-2.3.52/ncfsyn.m

   1 ## Copyright (C) 2011   Lukas F. Reichlin
   2 ##
   3 ## This file is part of LTI Syncope.
   4 ##
   5 ## LTI Syncope is free software: you can redistribute it and/or modify
   6 ## it under the terms of the GNU General Public License as published by
   7 ## the Free Software Foundation, either version 3 of the License, or
   8 ## (at your option) any later version.
   9 ##
  10 ## LTI Syncope is distributed in the hope that it will be useful,
  11 ## but WITHOUT ANY WARRANTY; without even the implied warranty of
  12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  13 ## GNU General Public License for more details.
  14 ##
  15 ## You should have received a copy of the GNU General Public License
  16 ## along with LTI Syncope.  If not, see <http://www.gnu.org/licenses/>.
  17
  18 ## -*- texinfo -*-
  19 ## @deftypefn{Function File} {[@var{K}, @var{N}, @var{gamma}, @var{info}] =} ncfsyn (@var{G}, @var{W1}, @var{W2}, @var{factor})
  20 ## Loop shaping H-infinity synthesis.  Compute positive feedback controller using
  21 ## the McFarlane/Glover normalized coprime factor (NCF) loop shaping design procedure.
  22 ##
  23 ## @strong{Inputs}
  24 ## @table @var
  25 ## @item G
  26 ## LTI model of plant.
  27 ## @item W1
  28 ## LTI model of precompensator.  Model must be SISO or of appropriate size.
  29 ## An identity matrix is taken if @var{W1} is not specified or if an empty model
  30 ## @code{[]} is passed.
  31 ## @item W2
  32 ## LTI model of postcompensator.  Model must be SISO or of appropriate size.
  33 ## An identity matrix is taken if @var{W2} is not specified or if an empty model
  34 ## @code{[]} is passed.
  35 ## @item factor
  36 ## @code{factor = 1} implies that an optimal controller is required.
  37 ## @code{factor > 1} implies that a suboptimal controller is required,
  38 ## achieving a performance that is @var{factor} times less than optimal.
  39 ## Default value is 1.
  40 ## @end table
  41 ##
  42 ## @strong{Outputs}
  43 ## @table @var
  44 ## @item K
  45 ## State-space model of the H-infinity loop-shaping controller.
  46 ## @item N
  47 ## State-space model of the closed loop depicted below.
  48 ## @item gamma
  49 ## L-infinity norm of @var{N}.  @code{gamma = norm (N, inf)}.
  50 ## @item info
  51 ## Structure containing additional information.
  52 ## @item info.emax
  53 ## Nugap robustness.  @code{emax = inv (gamma)}.
  54 ## @item info.Gs
  55 ## Shaped plant.  @code{Gs = W2 * G * W1}.
  56 ## @item info.Ks
  57 ## Controller for shaped plant.  @code{Ks = ncfsyn (Gs)}.
  58 ## @item info.rcond
  59 ## Estimates of the reciprocal condition numbers of the Riccati equations
  60 ## and a few other things.  For details, see the description of the
  61 ## corresponding SLICOT algorithm.
  62 ## @end table
  63 ##
  64 ## @strong{Block Diagram of N}
  65 ## @example
  66 ## @group
  67 ##
  68 ##             ^ z1              ^ z2
  69 ##             |                 |
  70 ##  w1  +      |   +--------+    |            +--------+
  71 ## ----->(+)---+-->|   Ks   |----+--->(+)---->|   Gs   |----+
  72 ##        ^ +      +--------+          ^      +--------+    |
  73 ##        |                        w2  |                    |
  74 ##        |                                                 |
  75 ##        +-------------------------------------------------+
  76 ## @end group
  77 ## @end example
  78 ##
  79 ## @strong{Algorithm}@*
  80 ## Uses SLICOT SB10ID, SB10KD and SB10ZD by courtesy of
  81 ## @uref{http://www.slicot.org, NICONET e.V.}
  82 ## @end deftypefn
  83
  84 ## Author: Lukas Reichlin <lukas.reichlin@gmail.com>
  85 ## Created: July 2011
  86 ## Version: 0.1
  87
  88 function [K, varargout] = ncfsyn (G, W1 = [], W2 = [], factor = 1.0)
  89
  90   if (nargin == 0 || nargin > 4)
  91     print_usage ();
  92   endif
  93
  94   if (! isa (G, "lti"))
  95     error ("ncfsyn: first argument must be an LTI system");
  96   endif
  97
  98   if (! is_real_scalar (factor) || factor < 1.0)
  99     error ("ncfsyn: fourth argument invalid");
 100   endif
 101
 102   [p, m] = size (G);
 103
 104   W1 = __adjust_weighting__ (W1, m);
 105   W2 = __adjust_weighting__ (W2, p);
 106
 107   Gs = W2 * G * W1;            # shaped plant
 108
 109   [a, b, c, d, tsam] = ssdata (Gs);
 110
 111   ## synthesis
 112   if (isct (Gs))               # continuous-time
 113     [ak, bk, ck, dk, rcond] = slsb10id (a, b, c, d, factor);
 114   elseif (any (d(:)))          # discrete-time, d != 0
 115     [ak, bk, ck, dk, rcond] = slsb10zd (a, b, c, d, factor, 0.0);
 116   else                         # discrete-time, d == 0
 117     [ak, bk, ck, dk, rcond] = slsb10kd (a, b, c, factor);
 118   endif
 119
 120   ## controller
 121   Ks = ss (ak, bk, ck, dk, tsam);
 122
 123   K = W1 * Ks * W2;
 124
 125   if (nargout > 1)
 126     ## FIXME: is this really the same thing as the dark side does?
 127     N = append (eye (p), Ks, Gs);
 128     M = [zeros(p,p), zeros(p,m),     eye(p);
 129              eye(p), zeros(p,m), zeros(p,p);
 130          zeros(m,p),     eye(m), zeros(m,p)];
 131     in_idx = [1:p, 2*p+(1:m)];
 132     out_idx = 1:p+m;
 133     N = mconnect (N, M, in_idx, out_idx);
 134     varargout{1} = N;
 135     if (nargout > 2)
 136       gamma = norm (N, inf);
 137       varargout{2} = gamma;
 138       if (nargout > 3)
 139         varargout{3} = struct ("emax", inv (gamma), "Gs", Gs, "Ks", Ks, "rcond", rcond);
 140       endif
 141     endif
 142   endif
 143
 144 endfunction
 145
 146
 147 function W = __adjust_weighting__ (W, s)
 148
 149   if (isempty (W))
 150     W = ss (eye (s));
 151   else
 152     W = ss (W);
 153     ## if (! isstable (W))
 154     ##   error ("ncfsyn: %s must be stable", inputname (1));
 155     ## endif
 156     ## if (! isminimumphase (W))
 157     ##   error ("ncfsyn: %s must be minimum-phase", inputname (1));
 158     ## endif
 159     [p, m] = size (W);
 160     if (m == s && p == s)      # model is of correct size
 161       return;
 162     elseif (m == 1 && p == 1)  # model is SISO
 163       tmp = W;
 164       for k = 2 : s
 165         W = append (W, tmp);   # stack SISO model s times
 166       endfor
 167     else                       # model is invalid
 168       error ("ncfsyn: %s must have 1 or %d inputs and outputs", inputname (1), s);
 169     endif
 170   endif
 171
 172 endfunction
 173
 174
 175 ## continuous-time case, direct access to sb10id
 176 %!shared AK, BK, CK, DK, RCOND, AKe, BKe, CKe, DKe, RCONDe
 177 %! A = [ -1.0  0.0  4.0  5.0 -3.0 -2.0
 178 %!       -2.0  4.0 -7.0 -2.0  0.0  3.0
 179 %!       -6.0  9.0 -5.0  0.0  2.0 -1.0
 180 %!       -8.0  4.0  7.0 -1.0 -3.0  0.0
 181 %!        2.0  5.0  8.0 -9.0  1.0 -4.0
 182 %!        3.0 -5.0  8.0  0.0  2.0 -6.0 ];
 183 %!
 184 %! B = [ -3.0 -4.0
 185 %!        2.0  0.0
 186 %!       -5.0 -7.0
 187 %!        4.0 -6.0
 188 %!       -3.0  9.0
 189 %!        1.0 -2.0 ];
 190 %!
 191 %! C = [  1.0 -1.0  2.0 -4.0  0.0 -3.0
 192 %!       -3.0  0.0  5.0 -1.0  1.0  1.0
 193 %!       -7.0  5.0  0.0 -8.0  2.0 -2.0 ];
 194 %!
 195 %! D = [  1.0 -2.0
 196 %!        0.0  4.0
 197 %!        5.0 -3.0 ];
 198 %!
 199 %! FACTOR = 1.0;
 200 %!
 201 %! [AK, BK, CK, DK, RCOND] = slsb10id (A, B, C, D, FACTOR);
 202 %!
 203 %! AKe = [ -39.0671    9.9293   22.2322  -27.4113   43.8655
 204 %!          -6.6117    3.0006   11.0878  -11.4130   15.4269
 205 %!          33.6805   -6.6934  -23.9953   14.1438  -33.4358
 206 %!         -32.3191    9.7316   25.4033  -24.0473   42.0517
 207 %!         -44.1655   18.7767   34.8873  -42.4369   50.8437 ];
 208 %!
 209 %! BKe = [ -10.2905  -16.5382  -10.9782
 210 %!          -4.3598   -8.7525   -5.1447
 211 %!           6.5962    1.8975    6.2316
 212 %!          -9.8770  -14.7041  -11.8778
 213 %!          -9.6726  -22.7309  -18.2692 ];
 214 %!
 215 %! CKe = [  -0.6647   -0.0599   -1.0376    0.5619    1.7297
 216 %!          -8.4202    3.9573    7.3094   -7.6283   10.6768 ];
 217 %!
 218 %! DKe = [  0.8466    0.4979   -0.6993
 219 %!         -1.2226   -4.8689   -4.5056 ];
 220 %!
 221 %! RCONDe = [ 0.13861D-01  0.90541D-02 ].';
 222 %!
 223 %!assert (AK, AKe, 1e-4);
 224 %!assert (BK, BKe, 1e-4);
 225 %!assert (CK, CKe, 1e-4);
 226 %!assert (DK, DKe, 1e-4);
 227 %!assert (RCOND, RCONDe, 1e-4);
 228
 229
 230 ## continuous-time case
 231 %!shared AK, BK, CK, DK, RCOND, AKe, BKe, CKe, DKe, RCONDe
 232 %! A = [ -1.0  0.0  4.0  5.0 -3.0 -2.0
 233 %!       -2.0  4.0 -7.0 -2.0  0.0  3.0
 234 %!       -6.0  9.0 -5.0  0.0  2.0 -1.0
 235 %!       -8.0  4.0  7.0 -1.0 -3.0  0.0
 236 %!        2.0  5.0  8.0 -9.0  1.0 -4.0
 237 %!        3.0 -5.0  8.0  0.0  2.0 -6.0 ];
 238 %!
 239 %! B = [ -3.0 -4.0
 240 %!        2.0  0.0
 241 %!       -5.0 -7.0
 242 %!        4.0 -6.0
 243 %!       -3.0  9.0
 244 %!        1.0 -2.0 ];
 245 %!
 246 %! C = [  1.0 -1.0  2.0 -4.0  0.0 -3.0
 247 %!       -3.0  0.0  5.0 -1.0  1.0  1.0
 248 %!       -7.0  5.0  0.0 -8.0  2.0 -2.0 ];
 249 %!
 250 %! D = [  1.0 -2.0
 251 %!        0.0  4.0
 252 %!        5.0 -3.0 ];
 253 %!
 254 %! FACTOR = 1.0;
 255 %!
 256 %! G = ss (A, B, C, D);
 257 %! K = ncfsyn (G, [], [], FACTOR);
 258 %! [AK, BK, CK, DK] = ssdata (K);
 259 %!
 260 %! AKe = [ -39.0671    9.9293   22.2322  -27.4113   43.8655
 261 %!          -6.6117    3.0006   11.0878  -11.4130   15.4269
 262 %!          33.6805   -6.6934  -23.9953   14.1438  -33.4358
 263 %!         -32.3191    9.7316   25.4033  -24.0473   42.0517
 264 %!         -44.1655   18.7767   34.8873  -42.4369   50.8437 ];
 265 %!
 266 %! BKe = [ -10.2905  -16.5382  -10.9782
 267 %!          -4.3598   -8.7525   -5.1447
 268 %!           6.5962    1.8975    6.2316
 269 %!          -9.8770  -14.7041  -11.8778
 270 %!          -9.6726  -22.7309  -18.2692 ];
 271 %!
 272 %! CKe = [  -0.6647   -0.0599   -1.0376    0.5619    1.7297
 273 %!          -8.4202    3.9573    7.3094   -7.6283   10.6768 ];
 274 %!
 275 %! DKe = [  0.8466    0.4979   -0.6993
 276 %!         -1.2226   -4.8689   -4.5056 ];
 277 %!
 278 %! RCONDe = [ 0.13861D-01  0.90541D-02 ];
 279 %!
 280 %!assert (AK, AKe, 1e-4);
 281 %!assert (BK, BKe, 1e-4);
 282 %!assert (CK, CKe, 1e-4);
 283 %!assert (DK, DKe, 1e-4);
 284
 285
 286 ## discrete-time case D==0, direct access to sb10kd
 287 %!shared AK, BK, CK, DK, RCOND, AKe, BKe, CKe, DKe, RCONDe
 288 %! A = [  0.2  0.0  0.3  0.0 -0.3 -0.1
 289 %!       -0.3  0.2 -0.4 -0.3  0.0  0.0
 290 %!       -0.1  0.1 -0.1  0.0  0.0 -0.3
 291 %!        0.1  0.0  0.0 -0.1 -0.1  0.0
 292 %!        0.0  0.3  0.6  0.2  0.1 -0.4
 293 %!        0.2 -0.4  0.0  0.0  0.2 -0.2 ];
 294 %!
 295 %! B = [ -1.0 -2.0
 296 %!        1.0  3.0
 297 %!       -3.0 -4.0
 298 %!        1.0 -2.0
 299 %!        0.0  1.0
 300 %!        1.0  5.0 ];
 301 %!
 302 %! C = [  1.0 -1.0  2.0 -2.0  0.0 -3.0
 303 %!       -3.0  0.0  1.0 -1.0  1.0 -1.0 ];
 304 %!
 305 %! FACTOR = 1.1;
 306 %!
 307 %! [AK, BK, CK, DK, RCOND] = slsb10kd (A, B, C, FACTOR);
 308 %!
 309 %! AKe = [  0.0337   0.0222   0.0858   0.1264  -0.1872   0.1547
 310 %!          0.4457   0.0668  -0.2255  -0.3204  -0.4548  -0.0691
 311 %!         -0.2419  -0.2506  -0.0982  -0.1321  -0.0130  -0.0838
 312 %!         -0.4402   0.3654  -0.0335  -0.2444   0.6366  -0.6469
 313 %!         -0.3623   0.3854   0.4162   0.4502   0.0065   0.1261
 314 %!         -0.0121  -0.4377   0.0604   0.2265  -0.3389   0.4542 ];
 315 %!
 316 %! BKe = [  0.0931  -0.0269
 317 %!         -0.0872   0.1599
 318 %!          0.0956  -0.1469
 319 %!         -0.1728   0.0129
 320 %!          0.2022  -0.1154
 321 %!          0.2419  -0.1737 ];
 322 %!
 323 %! CKe = [ -0.3677   0.2188   0.0403  -0.0854   0.3564  -0.3535
 324 %!          0.1624  -0.0708   0.0058   0.0606  -0.2163   0.1802 ];
 325 %!
 326 %! DKe = [ -0.0857  -0.0246
 327 %!          0.0460   0.0074 ];
 328 %!
 329 %! RCONDe = [ 0.11269D-01  0.17596D-01  0.18225D+00  0.75968D-03 ].';
 330 %!
 331 %!assert (AK, AKe, 1e-4);
 332 %!assert (BK, BKe, 1e-4);
 333 %!assert (CK, CKe, 1e-4);
 334 %!assert (DK, DKe, 1e-4);
 335 %!assert (RCOND, RCONDe, 1e-4);
 336
 337
 338 ## discrete-time case D==0
 339 %!shared AK, BK, CK, DK, RCOND, AKe, BKe, CKe, DKe, RCONDe
 340 %! A = [  0.2  0.0  0.3  0.0 -0.3 -0.1
 341 %!       -0.3  0.2 -0.4 -0.3  0.0  0.0
 342 %!       -0.1  0.1 -0.1  0.0  0.0 -0.3
 343 %!        0.1  0.0  0.0 -0.1 -0.1  0.0
 344 %!        0.0  0.3  0.6  0.2  0.1 -0.4
 345 %!        0.2 -0.4  0.0  0.0  0.2 -0.2 ];
 346 %!
 347 %! B = [ -1.0 -2.0
 348 %!        1.0  3.0
 349 %!       -3.0 -4.0
 350 %!        1.0 -2.0
 351 %!        0.0  1.0
 352 %!        1.0  5.0 ];
 353 %!
 354 %! C = [  1.0 -1.0  2.0 -2.0  0.0 -3.0
 355 %!       -3.0  0.0  1.0 -1.0  1.0 -1.0 ];
 356 %!
 357 %! FACTOR = 1.1;
 358 %!
 359 %! G = ss (A, B, C, [], 1);  # value of sampling time doesn't matter
 360 %! K = ncfsyn (G, [], [], FACTOR);
 361 %! [AK, BK, CK, DK] = ssdata (K);
 362 %!
 363 %! AKe = [  0.0337   0.0222   0.0858   0.1264  -0.1872   0.1547
 364 %!          0.4457   0.0668  -0.2255  -0.3204  -0.4548  -0.0691
 365 %!         -0.2419  -0.2506  -0.0982  -0.1321  -0.0130  -0.0838
 366 %!         -0.4402   0.3654  -0.0335  -0.2444   0.6366  -0.6469
 367 %!         -0.3623   0.3854   0.4162   0.4502   0.0065   0.1261
 368 %!         -0.0121  -0.4377   0.0604   0.2265  -0.3389   0.4542 ];
 369 %!
 370 %! BKe = [  0.0931  -0.0269
 371 %!         -0.0872   0.1599
 372 %!          0.0956  -0.1469
 373 %!         -0.1728   0.0129
 374 %!          0.2022  -0.1154
 375 %!          0.2419  -0.1737 ];
 376 %!
 377 %! CKe = [ -0.3677   0.2188   0.0403  -0.0854   0.3564  -0.3535
 378 %!          0.1624  -0.0708   0.0058   0.0606  -0.2163   0.1802 ];
 379 %!
 380 %! DKe = [ -0.0857  -0.0246
 381 %!          0.0460   0.0074 ];
 382 %!
 383 %! RCONDe = [ 0.11269D-01  0.17596D-01  0.18225D+00  0.75968D-03 ].';
 384 %!
 385 %!assert (AK, AKe, 1e-4);
 386 %!assert (BK, BKe, 1e-4);
 387 %!assert (CK, CKe, 1e-4);
 388 %!assert (DK, DKe, 1e-4);
 389
 390
 391 ## discrete-time case D!=0, direct access to sb10zd
 392 %!shared AK, BK, CK, DK, RCOND, AKe, BKe, CKe, DKe, RCONDe
 393 %! A = [  0.2  0.0  3.0  0.0 -0.3 -0.1
 394 %!       -3.0  0.2 -0.4 -0.3  0.0  0.0
 395 %!       -0.1  0.1 -1.0  0.0  0.0 -3.0
 396 %!        1.0  0.0  0.0 -1.0 -1.0  0.0
 397 %!        0.0  0.3  0.6  2.0  0.1 -0.4
 398 %!        0.2 -4.0  0.0  0.0  0.2 -2.0 ];
 399 %!
 400 %! B = [ -1.0 -2.0
 401 %!        1.0  3.0
 402 %!       -3.0 -4.0
 403 %!        1.0 -2.0
 404 %!        0.0  1.0
 405 %!        1.0  5.0 ];
 406 %!
 407 %! C = [  1.0 -1.0  2.0 -2.0  0.0 -3.0
 408 %!       -3.0  0.0  1.0 -1.0  1.0 -1.0
 409 %!        2.0  4.0 -3.0  0.0  5.0  1.0 ];
 410 %!
 411 %! D = [ 10.0 -6.0
 412 %!       -7.0  8.0
 413 %!        2.0 -4.0 ];
 414 %!
 415 %! FACTOR = 1.1;
 416 %!
 417 %! [AK, BK, CK, DK, RCOND] = slsb10zd (A, B, C, D, FACTOR, 0.0);
 418 %!
 419 %! AKe = [  1.0128   0.5101  -0.1546   1.1300   3.3759   0.4911
 420 %!         -2.1257  -1.4517  -0.4486   0.3493  -1.5506  -1.4296
 421 %!         -1.0930  -0.6026  -0.1344   0.2253  -1.5625  -0.6762
 422 %!          0.3207   0.1698   0.2376  -1.1781  -0.8705   0.2896
 423 %!          0.5017   0.9006   0.0668   2.3613   0.2049   0.3703
 424 %!          1.0787   0.6703   0.2783  -0.7213   0.4918   0.7435 ];
 425 %!
 426 %! BKe = [  0.4132   0.3112  -0.8077
 427 %!          0.2140   0.4253   0.1811
 428 %!         -0.0710   0.0807   0.3558
 429 %!         -0.0121  -0.2019   0.0249
 430 %!          0.1047   0.1399  -0.0457
 431 %!         -0.2542  -0.3472   0.0523 ];
 432 %!
 433 %! CKe = [ -0.0372  -0.0456  -0.0040   0.0962  -0.2059  -0.0571
 434 %!          0.1999   0.2994   0.1335  -0.0251  -0.3108   0.2048 ];
 435 %!
 436 %! DKe = [  0.0629  -0.0022   0.0363
 437 %!         -0.0228   0.0195   0.0600 ];
 438 %!
 439 %! RCONDe = [ 0.27949D-03  0.66679D-03  0.45677D-01  0.23433D-07  0.68495D-01  0.76854D-01 ].';
 440 %!
 441 %!assert (AK, AKe, 1e-4);
 442 %!assert (BK, BKe, 1e-4);
 443 %!assert (CK, CKe, 1e-4);
 444 %!assert (DK, DKe, 1e-4);
 445 %!assert (RCOND, RCONDe, 1e-4);
 446
 447
 448 ## discrete-time case D!=0
 449 %!shared AK, BK, CK, DK, RCOND, AKe, BKe, CKe, DKe, RCONDe
 450 %! A = [  0.2  0.0  3.0  0.0 -0.3 -0.1
 451 %!       -3.0  0.2 -0.4 -0.3  0.0  0.0
 452 %!       -0.1  0.1 -1.0  0.0  0.0 -3.0
 453 %!        1.0  0.0  0.0 -1.0 -1.0  0.0
 454 %!        0.0  0.3  0.6  2.0  0.1 -0.4
 455 %!        0.2 -4.0  0.0  0.0  0.2 -2.0 ];
 456 %!
 457 %! B = [ -1.0 -2.0
 458 %!        1.0  3.0
 459 %!       -3.0 -4.0
 460 %!        1.0 -2.0
 461 %!        0.0  1.0
 462 %!        1.0  5.0 ];
 463 %!
 464 %! C = [  1.0 -1.0  2.0 -2.0  0.0 -3.0
 465 %!       -3.0  0.0  1.0 -1.0  1.0 -1.0
 466 %!        2.0  4.0 -3.0  0.0  5.0  1.0 ];
 467 %!
 468 %! D = [ 10.0 -6.0
 469 %!       -7.0  8.0
 470 %!        2.0 -4.0 ];
 471 %!
 472 %! FACTOR = 1.1;
 473 %!
 474 %! G = ss (A, B, C, D, 1);  # value of sampling time doesn't matter
 475 %! K = ncfsyn (G, [], [], FACTOR);
 476 %! [AK, BK, CK, DK] = ssdata (K);
 477 %!
 478 %! AKe = [  1.0128   0.5101  -0.1546   1.1300   3.3759   0.4911
 479 %!         -2.1257  -1.4517  -0.4486   0.3493  -1.5506  -1.4296
 480 %!         -1.0930  -0.6026  -0.1344   0.2253  -1.5625  -0.6762
 481 %!          0.3207   0.1698   0.2376  -1.1781  -0.8705   0.2896
 482 %!          0.5017   0.9006   0.0668   2.3613   0.2049   0.3703
 483 %!          1.0787   0.6703   0.2783  -0.7213   0.4918   0.7435 ];
 484 %!
 485 %! BKe = [  0.4132   0.3112  -0.8077
 486 %!          0.2140   0.4253   0.1811
 487 %!         -0.0710   0.0807   0.3558
 488 %!         -0.0121  -0.2019   0.0249
 489 %!          0.1047   0.1399  -0.0457
 490 %!         -0.2542  -0.3472   0.0523 ];
 491 %!
 492 %! CKe = [ -0.0372  -0.0456  -0.0040   0.0962  -0.2059  -0.0571
 493 %!          0.1999   0.2994   0.1335  -0.0251  -0.3108   0.2048 ];
 494 %!
 495 %! DKe = [  0.0629  -0.0022   0.0363
 496 %!         -0.0228   0.0195   0.0600 ];
 497 %!
 498 %! RCONDe = [ 0.27949D-03  0.66679D-03  0.45677D-01  0.23433D-07  0.68495D-01  0.76854D-01 ].';
 499 %!
 500 %!assert (AK, AKe, 1e-4);
 501 %!assert (BK, BKe, 1e-4);
 502 %!assert (CK, CKe, 1e-4);
 503 %!assert (DK, DKe, 1e-4);