Artifact 4bd66676aa3250ce1fbad5d5c7ecf6802e7e6ac410a9c1c54ff6554b40793ffe:
- Executable file
mttroot/mtt/bin/trans/dae2cse_r
— part of check-in
[ea9834becc]
at
2002-06-05 11:14:51
on branch origin/optimise-algebraic-equations
— ae.r now generated using def2write_r like cse?.r
fix_c.r called at ese2rdae stage so that pow gets fixed in ae.r.These changes produce the desired result (optimised algebraic equations) but
have highlighted a problem; when optimisation fails, Reduce does not write
a result. For complicated systems, this can lead to missing assignments in
the resultant code. (user: geraint@users.sourceforge.net, size: 11101) [annotate] [blame] [check-ins using] [more...]
#! /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.20 2002/04/28 18:41:26 geraint ## Fixed [ 549658 ] awk should be gawk. ## Replaced calls to awk with call to gawk. ## ## Revision 1.19 2001/10/26 01:01:49 geraint ## fixcc when rdae_is_dae (-cr). ## ## Revision 1.18 2001/10/05 23:37:32 geraint ## Fixed assignment statement in ae.r when RHS=0. ## ## Revision 1.17 2001/07/27 23:29:10 geraint ## Optimises only when requested (-opt). ## ## Revision 1.16 2001/07/13 04:54:04 geraint ## Branch merge: numerical-algebraic-solution back to main. ## ## Revision 1.15.2.4 2001/06/26 00:55:48 geraint ## Writes algebraic equation Jacobian _aej.r (not used yet). ## ## Revision 1.15.2.3 2001/05/09 00:19:22 geraint ## Fixed EOF error when MTTNYZ=0. ## ## Revision 1.15.2.2 2001/05/05 20:50:16 geraint ## Fixed errors when MTTNx=0. ## ## Revision 1.15.2.1 2001/05/04 04:07:24 geraint ## Numerical solution of algebraic equations. ## sys_ae.cc written for unsolved inputs. ## Solution of equations using hybrd from MINPACK (as used by Octave fsolve). ## ## Revision 1.15 2001/03/19 02:28:52 geraint ## Branch merge: merging-ode2odes-exe back to MAIN. ## ## Revision 1.14.2.1 2001/03/19 00:29:08 geraint ## Parse switches (-A) before calling def2write_r. ## Update $1_def.* instead of removing. ## ## Revision 1.14 2000/12/28 12:24:35 peterg ## *** empty log message *** ## ## Revision 1.13 2000/10/11 08:52:46 peterg ## Creates csex (cse with dxe only) rep. ## ## Revision 1.12 2000/10/10 21:00:58 peterg ## New code genration ## ## Revision 1.11 1998/11/26 09:18:55 peterg ## Incluse subs.r ## ## 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 and optimise option solve=0; solve_msg='' optimise=''; optimise_msg='' while [ -n "`echo $1 | grep '^-'`" ]; do case $1 in -A ) solve=1 solve_msg=' with explicit solution of algebraic equations' ;; -fixcc ) fixcc=`echo 'in "'$MTT_LIB'/reduce/fix_c.r";'` fix_msg='fixing c and cc code'; ;; -optimise) optimise='-optimise' optimise_msg=' with optimisation' ;; *) echo "$1 is an invalid argument - ignoring" ;; esac shift done # Create the reduce output code def2write_r $optimise $1 ae def2write_r $optimise $1 cse def2write_r $optimise $1 csex # Version without E matrix def2write_r $optimise $1 cseo echo "Creating $1_ae.r $optimise_msg" echo "Creating $1_cse.r $solve_msg $optimise_msg $fix_msg" echo "Creating $1_csex.r $optimise_msg" echo "Creating $1_cseo.r $optimise_msg" # Remove the old log file rm -f dae2cse_r.log # Remove some files rm -f $1_ae.r? $1_cse.r? $1_cseo.r? # 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; %Fix c code if required $fixcc 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 dont 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 and $solve OUT "$1_aej.r"; IF (MTTNyz>0) THEN % as above BEGIN WRITE "MATRIX MTTyzj(",MTTNyz,",",MTTNyz,")"; WRITE "%File: $1_aej.r"; FOR i := 1:MTTNyz DO FOR j := 1:MTTNyz DO BEGIN didj := df(MTTyz(i,1),mkid('mttui,j)); IF (didj NEQ 0) THEN WRITE "MTTyzj(",i,",",j,") := ",didj," +0"; END; END; WRITE ";END;"; SHUT "$1_aej.r"; IF MTTNyz>0 THEN % not $solve or solution failed BEGIN OUT "$1_ae.r1"; write "MATRIX MTTYZ(", MTTNyz, ",", 1, ")$"; SHUT "$1_ae.r1"; OUT "$1_ae.r2"; write "%File: $1_ae.r"; in ("$1_ae_write.r"); write "END;"; SHUT "$1_ae.r2"; END; % Create the matrix declarations OUT "$1_cse.r1"; write "%"; IF (MTTNx > 0) THEN BEGIN write "MATRIX MTTEdx(", MTTNx, ",", 1, ")$"; write "MATRIX MTTE(", MTTNx, ",", MTTNx, ")$"; END; SHUT "$1_cse.r1"; OUT "$1_csex.r1"; write "%File:$1_csex.r1"; IF (MTTNx > 0) THEN write "MATRIX MTTEdx(", MTTNx, ",", 1, ")$"; SHUT "$1_csex.r1"; IF MTTNy>0 THEN BEGIN OUT "$1_cseo.r1"; write "MATRIX MTTY(", MTTNy, ",", 1, ")$"; SHUT "$1_cseo.r1"; END; %%Create the _cse.r file OUT "$1_cse.r2"; write "%File: $1_cse.r"; in ("$1_cse_write.r"); write "in ""$1_cseo.r"";"; write "END;"; % % 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.r2"; OUT "$1_csex.r2"; write "%File: $1_cse.r"; in ("$1_csex_write.r"); write "END;"; SHUT "$1_csex.r2"; %Write out the output equations IF MTTNy>0 THEN BEGIN OUT "$1_cseo.r2"; write "%File: $1_cseo.r"; in ("$1_cseo_write.r"); write "END;"; SHUT "$1_cseo.r2"; END; quit; EOF touch $1_ae.r1 $1_ae.r2 touch $1_cseo.r1 touch $1_cseo.r2 cat $1_ae.r1 $1_ae.r2 > $1_ae.r cat $1_cse.r1 $1_cse.r2 > $1_cse.r cat $1_csex.r1 $1_csex.r2 > $1_csex.r cat $1_cseo.r1 $1_cseo.r2 > $1_cseo.r if [ "$solve" = "1" ]; then echo "Setting MTTNyz=0 in $1_def.r and updating other $1_def files" gawk '{ 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 mv mtt_junk $1_def.r for file in `ls $1_def.*`; do if [ $file != $1_def.r ]; then ext=`echo $file|gawk -F\. '{print $2}'` mtt -q $1 def $ext fi done fi # Now invoke the standard error handling. mtt_error_r dae2cse_r.log