#! /bin/sh
######################################
##### Model Transformation Tools #####
######################################
# Bourne shell script: dae2cse_r
# Differential-algebraic equations to constrained-state equations
# P.J.Gawthrop 14 June 1991, 8 Aug 1991, 2 April 1992, 14 April 1994, 28 Dec 94
# Copyright (c) P.J.Gawthrop 1991, 1992, 1994.
###############################################################
## Version control history
###############################################################
## $Id$
## $Log$
## Revision 1.10 1998/11/18 13:50:29 peterg
## Removed writeing of EYz matrix
##
## Revision 1.9 1998/11/18 10:53:38 peterg
## Put in some more "IF MTTNx>0 THEN" to avoid error messages when no
## states.
##
## Revision 1.8 1998/11/10 08:54:34 peterg
## Put in "IF MTTNx>0 THEN" to prevent probs when Nx=0
## -- still a couple of apparent error messages - but answers now
## correct
##
## Revision 1.7 1998/10/05 10:46:15 peterg
## Commented out redundant MTTY := MTTY + MTTEyx*MTTEdX;
##
## Revision 1.6 1998/07/19 12:44:35 peterg
## Set MTTYz := 0 if the array is empty - avoids irritating error
## message.
##
## Revision 1.5 1998/05/20 15:23:26 peterg
## Put MTTYz := MTTYz outsise the BEGIN/END
##
## Revision 1.4 1998/05/20 15:13:09 peterg
## Writes out algebraic equations (if any).
##
## Revision 1.3 1998/03/03 09:02:46 peterg
## Replaced MTTEyx*MTTEdX + MTTEyu*MTTdu; term
##
## Revision 1.2 1997/08/26 08:22:36 peterg
## Changed
## MTTY := MTTY + MTTEyx*MTTdX + MTTEyu*MTTdu;
## to
## MTTY := MTTY + MTTEyx*MTTEdX + MTTEyu*MTTdu;
##
## This sorts out the problem when dz appears in the output equation.
##
## Revision 1.1 1997/08/26 08:20:18 peterg
## Initial revision
##
## Revision 1.2 1996/08/25 09:57:30 peter
## Sorted out bug when MTTNz=0
##
## Revision 1.1 1996/08/15 16:47:02 peter
## Initial revision
##
###############################################################
#Explicit solution option
solve=0
while [ -n "`echo $1 | grep '^-'`" ]; do
case $1 in
-A )
solve=1;;
*)
echo "$1 is an invalid argument - ignoring" ;;
esac
shift
done
if [ "$solve" = "1" ]; then
echo "Creating $1_cse.r (with explicit solution of algebraic equations)"
else
echo "Creating $1_cse.r"
fi
# Remove the old log file
rm -f dae2cse_r.log
# Use reduce to accomplish the transformation
$SYMBOLIC >dae2cse_r.log << EOF
%Read the formatting function
in "$MTTPATH/trans/reduce_matrix.r";
OFF Echo;
OFF Nat;
ON NERO;
in "$1_def.r";
MTTdxs := MTTdX; %Save the symbolic form of dX
in "$1_subs.r";
in "$1_dae.r";
%Create F_x, F_y matrices - assumming equations are
% linear in dZ
IF MTTNz>0 THEN
BEGIN
IF MTTNx>0 THEN
BEGIN
% Find MTTFx;
write "% Find MTTFx;";
matrix MTTFx(MTTNx,MTTNz);
FOR j := 1:MTTNz DO
BEGIN
dzj := MTTdZ(j,1);
FOR i := 1:MTTNx DO
MTTFx(i,j) := df(MTTdX(i,1), dzj, 1);
END;
END;
% Find MTTFy;
write "% Find MTTFy;";
matrix MTTFy(MTTNy,MTTNz);
FOR j := 1:MTTNz DO
BEGIN
dzj := MTTdZ(j,1);
FOR i := 1:MTTNy DO
MTTFy(i,j) := df(MTTy(i,1), dzj, 1);
END;
%Create G_x, G_u matrices
write "%Create G_x, G_u matrices ";
% Find MTTGx;
IF MTTNx>0 THEN
BEGIN
write "% Find MTTGx;";
matrix MTTGx(MTTNz,MTTNx);
FOR j := 1:MTTNx DO
BEGIN
xj := MTTX(j,1);
FOR i := 1:MTTNz DO
MTTGx(i,j) := df(MTTZ(i,1), xj, 1);
END;
END;
% Find MTTGu;
write "% Find MTTGu;";
matrix MTTGu(MTTNz,MTTNu);
FOR j := 1:MTTNu DO
BEGIN
uj := MTTu(j,1);
FOR i := 1:MTTNz DO
MTTGu(i,j) := df(MTTZ(i,1), uj, 1);
END;
%Create E matrices
write "%Create E matrices";
IF MTTNx>0 THEN
BEGIN
matrix MTTExx(MTTNx,MTTNx); MTTExx := MTTFx*MTTGx;
matrix MTTExu(MTTNx,MTTNu); MTTExu := MTTFx*MTTGu;
matrix MTTEyx(MTTNy,MTTNx); MTTEyx := MTTFy*MTTGx;
matrix MTTE(MTTNx,MTTNx); MTTE := MTTI - MTTExx;
END;
matrix MTTEyu(MTTNy,MTTNu); MTTEyu := MTTFy*MTTGu;
%% The following gets rid of the dZs; there must be a better way.
MTTdZ1 := 0;
MTTdZ2 := 0;
MTTdZ3 := 0;
MTTdZ4 := 0;
MTTdZ5 := 0;
MTTdZ6 := 0;
MTTdZ7 := 0;
MTTdZ8 := 0;
MTTdZ9 := 0;
MTTdZ10 := 0;
MTTdZ11 := 0;
MTTdZ12 := 0;
MTTdZ13 := 0;
MTTdZ14 := 0;
MTTdZ15 := 0;
MTTdZ16 := 0;
MTTdZ17 := 0;
MTTdZ18 := 0;
MTTdZ19 := 0;
IF MTTNx>0 THEN
BEGIN
MTTEdX := MTTdX; %Ie MTTEdX is MTTdX with the dz terms deleted ie EdX.
MTTdX := MTTdXs; %Restore the symbolic dX
%% Add on input derivative terms
MTTEdX := MTTEdX + MTTExu*MTTdu;
END;
%%%%%MTTY := MTTY + MTTEyx*MTTEdX;
%%% This causes the matrix mismatch
%%% MTTdXs and MTTdu need setting in _def.r file
MTTY := MTTY + MTTEyu*MTTdu;
IF MTTNx>0 THEN
MTTY := MTTY + MTTEyx*(MTTE^(-1))*MTTEdX;
END; %%of MTTNz>0
IF MTTNz=0 THEN
BEGIN
MTTEdX := MTTdX;
MTTE := MTTI;
END;
IF (MTTNyz>0) AND ($solve>0) THEN
BEGIN
%%%% Try and solve algebraic loops!!
%Create list of the relevant equations
MTT_eqns := {};
FOR i := 1:MTTNyz DO
MTT_eqns := append(MTT_eqns,{MTTyz(i,1)});
%Create list of the relevant unknowns
MTT_unknowns := {};
FOR i := 1:MTTNyz DO
MTT_unknowns := append(MTT_unknowns,{MTTUi(i,1)});
%Solve the algebraic equations symbolically
MTT_sol := solve(MTT_eqns,MTT_unknowns);
%The result seems to be in an extra list - I don't know why
% So remove the outer list with first.
% But only if more than one list element!
if MTTNyz>1 THEN
MTT_sol := first(MTT_sol);
%Substitute back into the equations
FOR i := 1:MTTNyz DO
BEGIN
MTT_sol_i := first(MTT_sol); MTT_sol := rest(MTT_sol);
set(lhs(MTT_sol_i),rhs(MTT_sol_i));
END;
% No algebraic variables left!
MTTNYz := 0;
END; % IF MTTNyz>0
%%Create the _cse.r file
OUT "$1_cse.r";
% State derivative
MTT_Matrix := MTTEdX$
MTT_Matrix_name := "MTTEdX"$
MTT_Matrix_n := MTTNx$
MTT_Matrix_m := 1$
Reduce_Matrix()$
% Output
MTT_Matrix := MTTY$
MTT_Matrix_name := "MTTY"$
MTT_Matrix_n := MTTNy$
MTT_Matrix_m := 1$
Reduce_Matrix()$
% Inputs
MTT_Matrix := MTTU$
MTT_Matrix_name := "MTTU"$
MTT_Matrix_n := MTTNu$
MTT_Matrix_m := 1$
Reduce_Matrix()$
MTT_Matrix := MTTdU$
MTT_Matrix_name := "MTTdU"$
MTT_Matrix_n := MTTNu$
MTT_Matrix_m := 1$
Reduce_Matrix()$
% E matrix
MTT_Matrix := MTTE$
MTT_Matrix_name := "MTTE"$
MTT_Matrix_n := MTTNx$
MTT_Matrix_m := MTTNx$
Reduce_Matrix()$
% Eyx matrix
MTT_Matrix := MTTEyx$
MTT_Matrix_name := "MTTEyx"$
MTT_Matrix_n := MTTNy$
MTT_Matrix_m := MTTNx$
%Reduce_Matrix()$
% Yz
MTT_Matrix := MTTYz$
MTT_Matrix_name := "MTTYz"$
MTT_Matrix_n := MTTNyz$
MTT_Matrix_m := 1$
Reduce_Matrix()$
write ";END;";
SHUT "$1_cse.r";
quit;
EOF
if [ "$solve" = "1" ]; then
echo "Setting MTTNyz=0 in $1_def.r and removing other $1_def files"
awk '{
if ($1=="MTTNyz")
print "MTTNyz := 0;"
else print $0
}' $1_def.r > mtt_junk
# Make sure it preserves the time stamp!!
# and remove dependent reps
touch -r $1_def.r mtt_junk
rm $1_def.*
mv mtt_junk $1_def.r
fi
# Now invoke the standard error handling.
mtt_error_r dae2cse_r.log