Skip to content

Commit

Permalink
improve dynamic state selection
Browse files Browse the repository at this point in the history 
 - fix initial state selection matrix A
 - improve LOG_DSS dumps
  • Loading branch information
@OpenModelica-Hudson
Willi Braun authored and OpenModelica-Hudson committed Nov 2, 2016
1 parent 70c7422 commit 786193c
Show file tree
Hide file tree
Showing 2 changed files with 72 additions and 50 deletions.
54 changes: 27 additions & 27 deletions Compiler/BackEnd/IndexReduction.mo
View file
Expand Up @@ -1596,7 +1596,7 @@ end addStateSets;
protected function generateStateSets
"author: Frenkel TUD 2013-01
generate the found state sets for the system"
input StateSets iTplLst "nStates,nStateCandidates,nUnassignedEquations,StateCandidates,ConstraintEqns,OtherVars,OtherEqns";
input StateSets iTplLst "level,nStates,nStateCandidates,nUnassignedEquations,StateCandidates,ConstraintEqns,OtherVars,OtherEqns";
input Integer iSetIndex;
input BackendDAE.Variables iVars;
input BackendDAE.EquationArray iEqns;
Expand Down Expand Up @@ -1630,21 +1630,21 @@ algorithm
b := intGt(rang,1);
// generate Set Vars
(_,crset,setVars,crA,aVars,tp,crJ,varJ) := getSetVars(oSetIndex,rang,nStateCandidates,nUnassignedEquations,level);
// add Equations
// set.x = set.A*set.statecandidates
// der(set.x) = set.A*der(set.candidates)
expcrstates := List.map(stateCandidates,BackendVariable.varExp);
expcrstatesstart := List.map(expcrstates,makeStartExp);
expcrdstates := List.map(expcrstates,makeder);
expcrset := List.map(crset,Expression.crefExp);
expcrdset := List.map(expcrset,makeder);
expcrA := Expression.crefExp(crA);
expcrA := DAE.CAST(tp,expcrA);
tyExpCrStates := DAE.T_ARRAY(DAE.T_REAL_DEFAULT,{DAE.DIM_INTEGER(nStateCandidates)},DAE.emptyTypeSource);
op := if b then DAE.MUL_MATRIX_PRODUCT(DAE.T_ARRAY(DAE.T_REAL_DEFAULT,{DAE.DIM_INTEGER(rang)}, DAE.emptyTypeSource)) else DAE.MUL_SCALAR_PRODUCT(DAE.T_REAL_DEFAULT);
mulAstates := DAE.BINARY(expcrA,op,DAE.ARRAY(tyExpCrStates,true,expcrstates));
(mulAstates,_) := Expression.extendArrExp(mulAstates,false);
mulAdstates := DAE.BINARY(expcrA,op,DAE.ARRAY(tyExpCrStates,true,expcrdstates));
// add Equations
// set.x = set.A*set.statecandidates
// der(set.x) = set.A*der(set.candidates)
expcrstates := List.map(stateCandidates,BackendVariable.varExp);
expcrstatesstart := List.map(expcrstates,makeStartExp);
expcrdstates := List.map(expcrstates,makeder);
expcrset := List.map(crset,Expression.crefExp);
expcrdset := List.map(expcrset,makeder);
expcrA := Expression.crefExp(crA);
expcrA := DAE.CAST(tp,expcrA);
tyExpCrStates := DAE.T_ARRAY(DAE.T_REAL_DEFAULT,{DAE.DIM_INTEGER(nStateCandidates)},DAE.emptyTypeSource);
op := if b then DAE.MUL_MATRIX_PRODUCT(DAE.T_ARRAY(DAE.T_REAL_DEFAULT,{DAE.DIM_INTEGER(rang)}, DAE.emptyTypeSource)) else DAE.MUL_SCALAR_PRODUCT(DAE.T_REAL_DEFAULT);
mulAstates := DAE.BINARY(expcrA,op,DAE.ARRAY(tyExpCrStates,true,expcrstates));
(mulAstates,_) := Expression.extendArrExp(mulAstates,false);
mulAdstates := DAE.BINARY(expcrA,op,DAE.ARRAY(tyExpCrStates,true,expcrdstates));
(mulAdstates,_) := Expression.extendArrExp(mulAdstates,false);
expset := if b then DAE.ARRAY(DAE.T_ARRAY(DAE.T_REAL_DEFAULT,{DAE.DIM_INTEGER(rang)},DAE.emptyTypeSource),true,expcrset) else listHead(expcrset);
expderset := if b then DAE.ARRAY(DAE.T_ARRAY(DAE.T_REAL_DEFAULT,{DAE.DIM_INTEGER(rang)},DAE.emptyTypeSource),true,expcrdset) else listHead(expcrdset);
Expand Down Expand Up @@ -4078,7 +4078,7 @@ protected function getSetVars
"author: Frenkel TUD 2012-12"
input Integer index;
input Integer setsize;
input Integer nStates;
input Integer nCandidates;
input Integer nCEqns;
input Integer level;
output DAE.ComponentRef crstates;
Expand All @@ -4100,16 +4100,16 @@ algorithm
oSetVars := BackendVariable.generateArrayVar(crstates,BackendDAE.STATE(1,NONE()),tp,NONE());
oSetVars := List.map1(oSetVars,BackendVariable.setVarFixed,false);
crset := List.map(oSetVars,BackendVariable.varCref);
tp := if intGt(setsize,1) then DAE.T_ARRAY(DAE.T_INTEGER_DEFAULT,{DAE.DIM_INTEGER(setsize),DAE.DIM_INTEGER(nStates)}, DAE.emptyTypeSource)
else DAE.T_ARRAY(DAE.T_INTEGER_DEFAULT,{DAE.DIM_INTEGER(nStates)}, DAE.emptyTypeSource);
realtp := if intGt(setsize,1) then DAE.T_ARRAY(DAE.T_REAL_DEFAULT,{DAE.DIM_INTEGER(setsize),DAE.DIM_INTEGER(nStates)}, DAE.emptyTypeSource)
else DAE.T_ARRAY(DAE.T_REAL_DEFAULT,{DAE.DIM_INTEGER(nStates)}, DAE.emptyTypeSource);
tp := if intGt(setsize,1) then DAE.T_ARRAY(DAE.T_INTEGER_DEFAULT,{DAE.DIM_INTEGER(setsize),DAE.DIM_INTEGER(nCandidates)}, DAE.emptyTypeSource)
else DAE.T_ARRAY(DAE.T_INTEGER_DEFAULT,{DAE.DIM_INTEGER(nCandidates)}, DAE.emptyTypeSource);
realtp := if intGt(setsize,1) then DAE.T_ARRAY(DAE.T_REAL_DEFAULT,{DAE.DIM_INTEGER(setsize),DAE.DIM_INTEGER(nCandidates)}, DAE.emptyTypeSource)
else DAE.T_ARRAY(DAE.T_REAL_DEFAULT,{DAE.DIM_INTEGER(nCandidates)}, DAE.emptyTypeSource);
ocrA := ComponentReference.joinCrefs(set,ComponentReference.makeCrefIdent("A",tp,{}));
oAVars := BackendVariable.generateArrayVar(ocrA,BackendDAE.VARIABLE(),tp,NONE());
oAVars := List.map1(oAVars,BackendVariable.setVarFixed,true);
// add start value A[i,j] = if i==j then 1 else 0 via initial equations
oAVars := List.map1(oAVars,BackendVariable.setVarStartValue,DAE.ICONST(0));
oAVars := setSetAStart(oAVars,1,1,setsize,{});
oAVars := setSetAStart(oAVars,1,1,nCandidates,{});
tp := if intGt(nCEqns,1) then DAE.T_ARRAY(DAE.T_REAL_DEFAULT,{DAE.DIM_INTEGER(nCEqns)}, DAE.emptyTypeSource) else DAE.T_REAL_DEFAULT;
ocrJ := ComponentReference.joinCrefs(set,ComponentReference.makeCrefIdent("J",tp,{}));
oJVars := BackendVariable.generateArrayVar(ocrJ,BackendDAE.VARIABLE(),tp,NONE());
Expand All @@ -4121,11 +4121,11 @@ protected function setSetAStart
input list<BackendDAE.Var> iVars;
input Integer n;
input Integer r;
input Integer nStates;
input Integer nCandidates;
input list<BackendDAE.Var> iAcc;
output list<BackendDAE.Var> oAcc;
algorithm
oAcc := match(iVars,n,r,nStates,iAcc)
oAcc := match(iVars,n,r,nCandidates,iAcc)
local
BackendDAE.Var v;
list<BackendDAE.Var> rest;
Expand All @@ -4135,10 +4135,10 @@ algorithm
equation
start = if intEq(n,r) then 1 else 0;
v = BackendVariable.setVarStartValue(v,DAE.ICONST(start));
n1 = if intEq(n,nStates) then 1 else (n+1);
r1 = if intEq(n,nStates) then (r+1) else r;
n1 = if intEq(n,nCandidates) then 1 else (n+1);
r1 = if intEq(n,nCandidates) then (r+1) else r;
then
setSetAStart(rest,n1,r1,nStates,v::iAcc);
setSetAStart(rest,n1,r1,nCandidates,v::iAcc);
end match;
end setSetAStart;

Expand Down
68 changes: 45 additions & 23 deletions SimulationRuntime/c/simulation/solver/stateset.c
View file
Expand Up @@ -77,7 +77,7 @@ void initializeStateSetPivoting(DATA *data)
unsigned int aid = 0;
modelica_integer *A = NULL;

/* go troug all state sets */
/* go trough all state sets */
for(i=0; i<data->modelData->nStateSets; i++)
{
set = &(data->simulationInfo->stateSetData[i]);
Expand All @@ -86,7 +86,7 @@ void initializeStateSetPivoting(DATA *data)

memset(A, 0, set->nCandidates*set->nStates*sizeof(modelica_integer));

/* initialize row and col indizes */
/* initialize row and col indices */
for(n=0; n<set->nDummyStates; n++)
set->rowPivot[n] = n;

Expand Down Expand Up @@ -191,9 +191,10 @@ static void getAnalyticalJacobianSet(DATA* data, threadData_t *threadData, unsig

if(ACTIVE_STREAM(LOG_DSS_JAC))
{
char *buffer = (char*)malloc(sizeof(char)*data->simulationInfo->analyticJacobians[jacIndex].sizeCols*10);
char *buffer = (char*)malloc(sizeof(char)*data->simulationInfo->analyticJacobians[jacIndex].sizeCols*20);

infoStreamPrint(LOG_DSS_JAC, 1, "jacobian %dx%d [id: %d]", data->simulationInfo->analyticJacobians[jacIndex].sizeRows, data->simulationInfo->analyticJacobians[jacIndex].sizeCols, jacIndex);

for(i=0; i<data->simulationInfo->analyticJacobians[jacIndex].sizeRows; i++)
{
buffer[0] = 0;
Expand Down Expand Up @@ -277,7 +278,7 @@ static int comparePivot(modelica_integer *oldPivot, modelica_integer *newPivot,
modelica_integer entry = (i < nDummyStates) ? 1: 2;
newEnable[ newPivot[i] ] = entry;
oldEnable[ oldPivot[i] ] = entry;
}
}

for(i=0; i<nCandidates; i++)
{
Expand All @@ -301,6 +302,44 @@ static int comparePivot(modelica_integer *oldPivot, modelica_integer *newPivot,
return ret;
}

/*! \fn printStateSelectionInfo
*
* function prints actually information about current state selection
*
* \param [in] [data]
* \param [in] [set]
*
* \author wbraun
*/
void printStateSelectionInfo(DATA *data, STATE_SET_DATA *set)
{
long k, l;

infoStreamPrint(LOG_DSS, 1, "Select %ld states from %ld candidates.", set->nStates, set->nCandidates);
for(k=0; k < set->nCandidates; k++)
{
infoStreamPrint(LOG_DSS, 0, "[%ld] cadidate %s", k+1, set->statescandidates[k]->name);
}
messageClose(LOG_DSS);

infoStreamPrint(LOG_DSS, 1, "Selected states");
{
unsigned int aid = set->A->id - data->modelData->integerVarsData[0].info.id;
modelica_integer *Adump = &(data->localData[0]->integerVars[aid]);
for(k=0; k < set->nStates; k++)
{
for(l=0; l < set->nCandidates; l++)
{
if (Adump[k*set->nCandidates+l] == 1)
{
infoStreamPrint(LOG_DSS, 0, "[%ld] %s", k+1, set->statescandidates[k]->name);
}
}
}
}
messageClose(LOG_DSS);
}

/*! \fn stateSelection
*
* function to select the actual states
Expand Down Expand Up @@ -329,29 +368,12 @@ int stateSelection(DATA *data, threadData_t *threadData, char reportError, int s
modelica_integer* oldColPivot = (modelica_integer*) malloc(set->nCandidates * sizeof(modelica_integer));
modelica_integer* oldRowPivot = (modelica_integer*) malloc(set->nDummyStates * sizeof(modelica_integer));


/* debug */
if(ACTIVE_STREAM(LOG_DSS))
{
infoStreamPrint(LOG_DSS, 1, "StateSelection Set %ld. Select %ld states from %ld candidates.", i, set->nStates, set->nCandidates);
for(k=0; k < set->nCandidates; k++)
{
infoStreamPrint(LOG_DSS, 0, "[%ld] cadidate %s", k+1, set->statescandidates[k]->name);
}
infoStreamPrint(LOG_DSS, 1, "StateSelection Set %ld at time = %f", i, data->localData[0]->timeValue);
printStateSelectionInfo(data, set);
messageClose(LOG_DSS);

infoStreamPrint(LOG_DSS, 0, "StateSelection Matrix A");
{
unsigned int aid = set->A->id - data->modelData->integerVarsData[0].info.id;
modelica_integer *Adump = &(data->localData[0]->integerVars[aid]);
for(k=0; k < set->nCandidates; k++)
{
for(l=0; l < set->nStates; l++)
{
infoStreamPrint(LOG_DSS, 0, "A[%ld,%ld] = %ld", k+1, l+1, Adump[k*set->nCandidates+l]);
}
}
}
}
/* generate jacobian, stored in set->J */
getAnalyticalJacobianSet(data, threadData, i);
Expand Down

0 comments on commit 786193c

Please sign in to comment.