# Copyright (C) 1998 Kai P. Mueller
#
# This file is part of Octave. 
#
# Octave 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 2, or (at your option) any 
# later version. 
# 
# Octave 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 Octave; see the file COPYING.  If not, write to the Free 
# Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. 
 
function [sys] = buildssic(Clst,Ulst,Olst,Ilst,s1,s2,s3,s4,s5,s6,s7,s8)
#
# [sys] = buildssic(Clst,Ulst,Olst,Ilst,s1,s2,s3,s4,s5,s6,s7,s8)
#
# Form an arbitrary complex (open or closed loop) system in
# state-space form from several systems. "buildssic" can
# easily (despite it's cryptic syntax) integrate transfer functions
# from a complex block diagram into a single system with one call.
# This function is especially useful for building open loop
# interconnections for H_infinity and H2 designs or for closing
# loops with these controllers.
#
# Although this function is general purpose, the use of "sysgroup"
# "sysmult", "sysconnect" and the like ist recommended for standard
# operations since they can handle mixed discrete and continuous
# systems and also the names of inputs, outputs, and states.
# 
# The parameters consist of 4 lists which describe the connections
# outputs and inputs and up to 8 systems s1-s8.
# Format of the lists:
#
#     Clst: connection list, describes the input signal of
#           each system. The maximum number of rows of Clst is
#           equal to the sum of all inputs of s1-s8.
#           Example:
#             [1 2 -1; 2 1 0] ==> new input 1 is old inpout 1
#             + output 2 - output 1, new input 2 is old input 2
#             + output 1. The order of rows is arbitrary.
#
#     Ulst: if not empty the old inputs in vector Ulst will
#           be appended to the outputs. You need this if you
#           want to "pull out" the input of a system. Elements
#           are input numbers of s1-s8.
#
#     Olst: output list, specifiy the outputs of the resulting
#           systems. Elements are output numbers of s1-s8.
#           The numbers are alowed to be negative and may
#           appear in any order. An empty matrix means
#           all outputs.
#
#     Ilst: input list, specifiy the inputs of the resulting
#           systems. Elements are input numbers of s1-s8.
#           The numbers are alowed to be negative and may
#           appear in any order. An empty matrix means
#           all inputs.
#
# Example:  Very simple closed loop system.
#
#       w        e  +-----+   u  +-----+
#        --->o--*-->|  K  |--*-->|  G  |--*---> y
#            ^  |   +-----+  |   +-----+  |
#          - |  |            |            |
#            |  |            +----------------> u
#            |  |                         |
#            |  +-------------------------|---> e
#            |                            |
#            +----------------------------+
#
# The closed loop system GW can be optained by
#
#     GW = buildssic([1 2; 2 -1], [2], [1 2 3], [2], G, K);
#
# Clst: (1. row) connect input 1 (G) with output 2 (K).
#       (2. row) connect input 2 (K) with neg. output 1 (G).
# Ulst: append input of (2) K to the number of outputs.
# Olst: Outputs are output of 1 (G), 2 (K) and appended
#       output 3 (from Ulst).
# Ilst: the only input is 2 (K).
#
# Here is a real example:
#	                         +----+
#	    -------------------->| W1 |---> v1
#	z   |                    +----+
#	----|-------------+                   || GW   ||     => min.
#	    |             |                        vz   infty
#	    |    +---+    v      +----+
#	    *--->| G |--->O--*-->| W2 |---> v2
#	    |    +---+       |   +----+
#	    |                |
#           |                v
#          u                  y
#
# The closed loop system GW from [z; u]' to [v1; v2; y]' can be
# obtained by (all SISO systems):
#
#     GW = buildssic([1 4;2 4;3 1],[3],[2 3 5],[3 4],G,W1,W2,One);
#
# where "One" is a unity gain (auxillary) function with order 0.
# (e.g. One = ugain(1);)
#

# Written by Kai Mueller April 1998

  if((nargin < 5) || (nargin > 12))
    usage("[sys] = buildssic(Clst,Ulst,Olst,Ilst,s1,s2,s3,s4,s5,s6,s7,s8)");
  endif
  if (nargin >= 5)
    if (!is_struct(s1))
      error("---> s1 must be a structed system.");
    endif
    s1 = sysupdate(s1, "ss");
    [n, nz, m, p] = sysdimensions(s1);
    if (!n && !nz)
      error("---> pure static system must not be the first in list.");
    endif
    if (n && nz)
      error("---> cannot handle mixed continuous and discrete systems.");
    endif
    D_SYS = (nz > 0);
    [A,B,C,D,tsam] = sys2ss(s1);
    nt = n + nz;
  endif
  for ii = 6:nargin
    eval(["ss = s", num2str(ii-4), ";"]);
    if (!is_struct(ss))
      error("---> Parameter must be a structed system.");
    endif
    ss = sysupdate(ss, "ss");
    [n1, nz1, m1, p1] = sysdimensions(ss);
    if (n1 && nz1)
      error("---> cannot handle mixed continuous and discrete systems.");
    endif
    if (D_SYS)
      if (n1)
      	error("---> cannot handle mixed cont. and discr. systems.");
      endif
      if (tsam != sysgettsam(ss))
	error("---> sampling time of all systems must match.");
      endif
    endif
    [as,bs,cs,ds] = sys2ss(ss);
    nt1 = n1 + nz1;
    if (!nt1)
      # pure gain (pad B, C with zeros)
      B = [B, zeros(nt,m1)];
      C = [C; zeros(p1,nt)];
    else
      A = [A, zeros(nt,nt1); zeros(nt1,nt), as];
      B = [B, zeros(nt,m1);  zeros(nt1,m),  bs];
      C = [C, zeros(p,nt1);  zeros(p1,nt),  cs];
    endif
    D = [D, zeros(p,m1); zeros(p1,m), ds];
    n = n + n1;
    nz = nz + nz1;
    nt = nt + nt1;
    m = m + m1;
    p = p + p1;
  endfor

  # check maximum dimensions
  [nx, mx] = size(Clst);
  if (nx > m)
    error("---> more rows in Clst than total number of inputs.");
  endif
  if (mx > p+1)
    error("---> more cols in Clst than total number of outputs.");
  endif
  # empty vector Ulst is OK
  lul = length(Ulst);
  if (lul)
    if (!is_vector(Ulst))
      error("---> Input u list Ulst must be a vector.");
    endif
    if (lul > m)
      error("---> more values in Ulst than number of inputs.");
    endif
  endif
  if (!length(Olst))  Olst = [1:(p+lul)];  endif
  if (!length(Ilst))  Ilst = [1:m];        endif
  if (!is_vector(Olst))
    error("---> Output list Olst must be a vector.");
  endif
  if (!is_vector(Ilst))
    error("---> Input list Ilst must be a vector.");
  endif

  # build the feedback "K" from the interconnection data Clst
  K = zeros(m, p);
  inp_used = zeros(m,1);
  for ii = 1:nx
    xx = Clst(ii,:);
    iu = xx(1);
    if ((iu < 1) || (iu > m))
      error("---> Illegal value in first col of Clst.");
    endif
    if (inp_used(iu))
      error("---> Input specified more than once.");
    endif
    inp_used(iu) = 1;
    for kk = 2:mx
      it = xx(kk);
      if (abs(it) > p)
      	error("---> Illegal row value in Clst.");
      elseif (it)
	K(iu,abs(it)) = sign(it);
      endif
    endfor
  endfor

  # form the "closed loop", i.e replace u in
  # .
  # x = Ax + Bu
  #                            ~
  # y = Cx + Du   by   u = K*y+u
  #
  #            -1
  # R = (I-D*K)   must exist.
  # 
  R = eye(p) - D*K;
  if (rank(R) < m)
    error("---> singularity in algebraic loop.");
  else
    R = inv(R);
  endif
  A = A + B*K*R*C;
  B = B + B*K*R*D;
  C = R*C;
  D = R*D;

  # append old inputs u to the outputs (if lul > 0)
  kc = K*C;
  kdi = eye(m) + K*D;
  for ii = 1:lul
    it = Ulst(ii);
    if ((it < 1) || (it > m))
      error("---> Illegal value in Ulst.");
    endif
    C = [C; kc(it,:)];
    D = [D; kdi(it,:)];
  endfor

  # select and rearrange outputs
  nn = length(A);
  lol = length(Olst);
  Cnew = zeros(lol,nn);
  Dnew = zeros(lol,m);
  for ii = 1:lol
    iu = Olst(ii);
    if (!iu || (abs(iu) > p+lul))
      error("---> Illegal value in Olst.");
    endif
    Cnew(ii,:) = sign(iu)*C(abs(iu),:);
    Dnew(ii,:) = sign(iu)*D(abs(iu),:);
  endfor
  C = Cnew;
  D = Dnew;
  lil = length(Ilst);
  Bnew = zeros(nn,lil);
  Dnew = zeros(lol,lil);
  for ii = 1:lil
    iu = Ilst(ii);
    if (!iu || (abs(iu) > m))
      error("---> Illegal value in Ilst.");
    endif
    Bnew(:,ii) = sign(iu)*B(:,abs(iu));
    Dnew(:,ii) = sign(iu)*D(:,abs(iu));
  endfor

  sys = ss2sys(A, Bnew, C, Dnew, tsam, n, nz);

endfunction