|
|
|
|
||
|
DEFINITION MODULE MonWeathGen; (******************************************************************* Module MonWeathGen (Version 1.0) Copyright (c) 1996-2006 by Dimitrios Gyalistras and ETH Zurich. Purpose Stochastic simulation of monthly weather variables. Remarks Used is the multivariate AR1-model: X(k+1) = A(p)*X(k) + B(p)*E(k) Y(k+1) = Z(k+1)*S'(p) + M(p) where +Exp[ Fi(p)*Xi(k) + Gi(p) ] + Hi(p) if SKEW[ Yi(p) ] > 0 Zi(k) = { Xi(k) if SKEW[ Yi(p) ] small -Exp[ Fi(p)*Xi(k) + Gi(p) ] + Hi(p) if SKEW[ Yi(p) ] < 0 and k : time (time step = 1 month) p : phase within the year, i.e. p = (k+12) MOD 12 i : i-th element of a vector ' : denotes transposition of a vector X : state vector (dimension N) A : system matrix (NxN) B : input matrix (NxN, symetric) E : input vector of N independent random components from a N-dimensional normal distribution N(0,1). Y : output vector M : vector of expected values of Y S : vector of standard deviations of Y Z : vector of auxiliary variables with normal or log-normal distributions, i.e. (-Gi + Log [ +Zi - Hi ])/Fi if SKEW[ Yi ] >0 Xi = { Zi if SKEW[ Yi ] small (-Gi + Log [ -Zi - Hi ])/Fi if SKEW[ Yi ] <0 where Xi is normally distributed and Fi, Gi, and Hi are the parameters of the log-normally distributed variable Zi. Note: SKEW[*] = E[ (x-MU[*])^3 ]/SIG[*]^3 denotes the skewness of the random variable *. See also modules CycloStat and LogNTransf. Programming o Design Dimitrios Gyalistras 07/12/1996 o Implementation Dimitrios Gyalistras 07/12/1996 ETH Zurich Systems Ecology CHN E 35.1 Universitaetstrasse 16 8092 Zurich SWITZERLAND URLs: <mailto:RAMSES@env.ethz.ch> <http://www.sysecol.ethz.ch> <http://www.sysecol.ethz.ch/SimSoftware/RAMSES> Last revision of definition: 24/05/2002 AF *******************************************************************) IMPORT Errors; FROM LgMatrices IMPORT LMatrix; FROM LogNTransf IMPORT LnTrfType; CONST (* result codes returned *) allOk = Errors.allOk; notDone = Errors.onlyAnInsert; CONST MaxMWGVars = 32; Jan = 1; Dec = 12; TYPE Month = [Jan..Dec]; TYPE MWGLRealVec = ARRAY [1..MaxMWGVars] OF LONGREAL; MWGIntVec = ARRAY [1..MaxMWGVars] OF INTEGER; MWGLnTrfTVec = ARRAY [1..MaxMWGVars] OF LnTrfType; MWGLRealVecArr = ARRAY [Jan..Dec] OF MWGLRealVec; MWGLnTrfTVecArr = ARRAY [Jan..Dec] OF MWGLnTrfTVec; MWGLMatrixArr = ARRAY [Jan..Dec] OF LMatrix; (*****************************************************************************) (**************************************************) (*##### 1 - Allocation of needed objects #####*) (**************************************************) PROCEDURE ResetMWGLMatrixArr( VAR mwgLma: MWGLMatrixArr ); (* Sets all elements of mwgLma to notAllocatedLMatrix *) PROCEDURE AllocMWGLMatrixArr( VAR mwgLma : MWGLMatrixArr; nVars : INTEGER; VAR resCode : INTEGER; VAR errTxt : ARRAY OF CHAR ); (* Allocates all mwgLma matrices with size nVars x nVars and initializes them with DMConversions.UndefLONGREAL *) PROCEDURE DeallocMWGLMatrixArr( VAR mwgLma: MWGLMatrixArr ); (* Deallocates all matrices in mwgLma *) (***************************************************************************) (*##### 2 - Get/Set all parameters defining stochastic simulation #####*) (***************************************************************************) PROCEDURE MWGParamsHaveBeenSet():BOOLEAN; (* Returns true if SetMWGParams has been called succesfully at least once. *) PROCEDURE GetMWGParams( VAR nVars : INTEGER; VAR mu : MWGLRealVecArr; (* M(p) *) VAR sig : MWGLRealVecArr; (* S(p) *) VAR lnTrfType : MWGLnTrfTVecArr; (* noLnTrf, pos/negSkewLnTrf *) VAR lnTrfMu : MWGLRealVecArr; (* G(p) *) VAR lnTrfSig : MWGLRealVecArr; (* F(p) *) VAR lnTrfTheta : MWGLRealVecArr; (* H(p) *) VAR A, B : MWGLMatrixArr; (* A, B *) VAR allowedMin : MWGLRealVecArr; (* lower range for Yi(p) *) VAR allowedMax : MWGLRealVecArr; (* upper range for Yi(p) *) VAR maxIters : INTEGER; VAR resCode : INTEGER; VAR errTxt : ARRAY OF CHAR ); PROCEDURE SetMWGParams( nVars : INTEGER; VAR mu : MWGLRealVecArr; (* VAR for speed-up only *) VAR sig : MWGLRealVecArr; (* VAR for speed-up only *) VAR lnTrfType : MWGLnTrfTVecArr; (* VAR for speed-up only *) VAR lnTrfMu : MWGLRealVecArr; (* VAR for speed-up only *) VAR lnTrfSig : MWGLRealVecArr; (* VAR for speed-up only *) VAR lnTrfTheta : MWGLRealVecArr; (* VAR for speed-up only *) VAR A, B : MWGLMatrixArr; (* VAR for speed-up only *) VAR allowedMin : MWGLRealVecArr; (* VAR for speed-up only *) VAR allowedMax : MWGLRealVecArr; (* VAR for speed-up only *) maxIters : INTEGER; VAR resCode : INTEGER; VAR errTxt : ARRAY OF CHAR ); (* If during simulation an Yi(k+1) is found outside the range [allowedMin(p,i)..allowedMax(p,i)] and maxIters<>0, then the calculation of Y(k+1) is repeated up to maxIters times until all Yi are within their respective range. If maxIters have been reached variable wgOutputOK (see below) will return FALSE. If maxIters=0 all elements of Y(k+1) are adjusted to stay within their respective range (brute force method) and mwgOutputOK will always return TRUE. *) (**************************************************) (*##### 3 - Get/Set selected parameters. #####*) (**************************************************) (* NOTE: Calling one of the following procedures has an effect only if SetMWGParams has been called at least once before. The procedures have been optimized for efficiency, i.e. assignments are attempted without any further checks. In particular it is the client's responsibility to pass the correct MWGLMatrixArr-variables. *) PROCEDURE NumMWGVars():INTEGER; PROCEDURE GetMWGMu ( VAR mu : MWGLRealVecArr ); PROCEDURE SetMWGMu ( VAR mu : MWGLRealVecArr ); (* VAR for speed-up only *) PROCEDURE GetMWGSig ( VAR sig : MWGLRealVecArr ); PROCEDURE SetMWGSig ( VAR sig : MWGLRealVecArr ); (* VAR for speed-up only *) PROCEDURE GetMWGLnTrfPars( VAR lnTrfMu : MWGLRealVecArr; VAR lnTrfSig : MWGLRealVecArr; VAR lnTrfTheta : MWGLRealVecArr ); PROCEDURE SetMWGLnTrfPars( VAR lnTrfMu : MWGLRealVecArr; (* VAR for speed-up only *) VAR lnTrfSig : MWGLRealVecArr; (* VAR for speed-up only *) VAR lnTrfTheta : MWGLRealVecArr ); (* VAR for speed-up only *) PROCEDURE GetMWGMatrices ( VAR A, B : MWGLMatrixArr ); PROCEDURE SetMWGMatrices ( VAR A, B : MWGLMatrixArr ); (* VAR for speed-up only *) (******************************************************) (*##### 4 - Control of stochastic simulation #####*) (******************************************************) PROCEDURE InitMWG( VAR yInit : MWGLRealVec; (* initial state, VAR for speed-up only *) month : Month; (* month to which initial state refers *) initRandGen : BOOLEAN; (* if true, random number generator will be initialized *) z0,z1,z2 : INTEGER; (* seeds used to initialize random number generator *) VAR resCode : INTEGER; VAR errTxt : ARRAY OF CHAR ); PROCEDURE GetLastMWGOutput( VAR yOld: MWGLRealVec ); PROCEDURE ClcNewMWGOutput ( VAR yNew: MWGLRealVec ); VAR mwgOutputOK: BOOLEAN; (* Returns true if last call to ClcNewMWGOutput was succesful. Read only. *) END MonWeathGen.
|
||
|
|
|