- added some bugfixes for initialization and update some testcase to…

… currect values

git-svn-id: https://openmodelica.org/svn/OpenModelica/trunk@8667 f25d12d1-65f4-0310-ae8a-bbce733d8d8e

Compiler/BackEnd/BackendDAEUtil.mo

@@ -435,7 +435,7 @@ algorithm
         ((nie2,_)) = BackendEquation.traverseBackendDAEEqns(removedEqs,countInitialEqns,(nie1,whenClauseLst));
         dae1 = checkInitialSystem1(unfixed1,nie2,dae,funcs,vars,varsws,states,statesws);
       then
-        dae;
+        dae1;
 
     case (dae,_)
       equation

Compiler/susan_codegen/SimCode/SimCodeC.tpl

@@ -800,6 +800,7 @@ template functionInitializeDataStruc()
     }
 
     returnData->initFixed = init_fixed;
+    returnData->var_attr = var_attr;
     returnData->modelName = model_name;
     returnData->modelFilePrefix = model_fileprefix;
     returnData->statesNames = state_names;
@@ -1130,7 +1131,7 @@ template functionInitialResidual(list<SimEqSystem> residualEquations)
         let expPart = daeExp(exp, contextOther, &preExp /*BUFC*/,
                            &varDecls /*BUFD*/)
         '<%preExp%>localData->initialResiduals[i++] = <%expPart%>;
-if (sim_verbose >=LOG_INIT) { printf(" Residual[%d] : <%ExpressionDump.printExpStr(exp)%> = %f\n",i,localData->initialResiduals[i-1]); }'
+if (sim_verbose == LOG_RES_INIT) { printf(" Residual[%d] : <%ExpressionDump.printExpStr(exp)%> = %f\n",i,localData->initialResiduals[i-1]); }'
     ;separator="\n")
   <<
   int initial_residual()

c_runtime/simulation_init.cpp

@@ -30,6 +30,7 @@
 
 #include "simulation_init.h"
 #include "simulation_runtime.h"
+#include "solver_main.h"
 #include <math.h>
 
 /*
@@ -42,26 +43,18 @@ void leastSquare(long *nz, double *z, double *funcValue)
   int startIndPar = 2*globalData->nStates+globalData->nAlgebraic+globalData->intVariables.nAlgebraic+globalData->boolVariables.nAlgebraic;
 
  for (ind=0, indAct=0, indz=0; ind<globalData->nStates; ind++)
-    if (globalData->initFixed[indAct++]==0)
+    if (globalData->initFixed[indAct++]==0 )
           globalData->states[ind] = z[indz++];
 
   // for real parameters 
-  for (ind=0,indAct=startIndPar; ind<globalData->nParameters; ind++)
-    if (globalData->initFixed[indAct++]==0)
+  for (ind=0,indAct=startIndPar; ind<globalData->nParameters; ind++, indAct++)
+    if (globalData->initFixed[indAct]==0 and globalData->var_attr[indAct-globalData->nStates]==1)
       globalData->parameters[ind] = z[indz++];
 
   bound_parameters();
   functionODE();
   functionAlgebraics();
 
-/*  for (ind=0,indy=0,indAct=2*globalData->nStates; ind<globalData->nAlgebraic; ind++)
-    if (globalData->initFixed[indAct++]==1)
-      globalData->algebraics [ind] = static_y[indy++];
-
-      Comment from Bernhard: Even though algebraic variables are "fixed", they are calculated from
-      the states, so they should be allowed to change when states vary,
-      and NOT be replaced by their initial values as above.
-*/
   initial_residual();
 
   for (ind=0, *funcValue=0; ind<globalData->nInitialResiduals; ind++)
@@ -88,7 +81,7 @@ int reportResidualValue(double funcValue)
         cout << "residual[" << i << "] = " << globalData->initialResiduals[i] << endl;
       }
     }
-    return 0 /*-1*/;
+    return 1;
   }
   return 0;
 }
@@ -136,23 +129,23 @@ int simplex_initialization(long& nz,double *z)
   /* Start with stepping .5 in each direction. */
   for (ind=0;ind<nz;ind++)
   {
-    STEP[ind] = 0.5;
+    STEP[ind] = 1.0;
     VAR[ind]  = 0.0;
   }
 
    double STOPCR,SIMP;
    long IPRINT, NLOOP,IQUAD,IFAULT,MAXF;
 //C  Set max. no. of function evaluations = 5000, print every 100.
 
-      MAXF = 20 * nz;
+      MAXF = 50 * nz;
       IPRINT = sim_verbose >= LOG_INIT ? 100 : -1;
 
 //C  Set value for stopping criterion.   Stopping occurs when the
 //C  standard deviation of the values of the objective function at
 //C  the points of the current simplex < stopcr.
 
       STOPCR = 1.e-3;
-      NLOOP = 2*MAXF;//2*nz;
+      NLOOP = 2*MAXF;
 
 //C  Fit a quadratic surface to be sure a minimum has been found.
 
@@ -203,23 +196,15 @@ int simplex_initialization(long& nz,double *z)
  * residual of all equations (both continuous time eqns and initial eqns).
  */
 
-int initialize(const std::string*method)
+int initialize(const std::string init_method)
 {
   long nz;
   int ind, indAct, indz;
-  std::string init_method;
-
-  if (method == NULL) {
-   // init_method = std::string("newuoa");
-    init_method = std::string("simplex");
-  } else {
-    init_method = *method;
-  }
 
   for (ind=0, nz=0; ind<globalData->nStates; ind++){
     if (globalData->initFixed[ind]==0){
         if (sim_verbose >= LOG_INIT)
-            printf("Variable %s is unfixed.\n",globalData->statesNames[ind].name);
+            printf("State %s is unfixed.\n",globalData->statesNames[ind].name);
 
         nz++;
     }
@@ -228,9 +213,9 @@ int initialize(const std::string*method)
   int startIndPar = 2*globalData->nStates+globalData->nAlgebraic+globalData->intVariables.nAlgebraic+globalData->boolVariables.nAlgebraic;
   int endIndPar = startIndPar+globalData->nParameters;
   for (ind=startIndPar;ind<endIndPar; ind++){
-    if (globalData->initFixed[ind]==0){
+    if (globalData->initFixed[ind]==0 and globalData->var_attr[ind-globalData->nStates]==1){
         if (sim_verbose >= LOG_INIT)
-          printf("Variable %s is unfixed.\n",globalData->parametersNames[ind-startIndPar].name);
+          printf("Parameter %s is unfixed.\n",globalData->parametersNames[ind-startIndPar].name);
         nz++;
     }
   }
@@ -262,7 +247,7 @@ int initialize(const std::string*method)
   }
   // for real parameters
   for (ind=0,indAct=startIndPar; ind<globalData->nParameters; ind++) {
-    if (globalData->initFixed[indAct++]==0)
+    if (globalData->initFixed[indAct++]==0 and globalData->var_attr[indAct-globalData->nStates]==1)
       z[indz++] =  globalData->parameters[ind];
   }
 
@@ -280,3 +265,89 @@ int initialize(const std::string*method)
   return retVal;
 }
 
+int
+main_initialize(const std::string* method)
+{
+
+  std::string init_method;
+
+  if (method == NULL)
+    {
+      init_method = std::string("simplex");
+    }
+  else
+    {
+      init_method = *method;
+    }
+
+  saveall();
+  initial_function();
+  update_DAEsystem();
+  if (sim_verbose >= LOG_SOLVER)
+    {
+      sim_result->emit();
+    }
+  // do some initial stuff
+  globalData->init = 1;
+  initial_function();
+  if (sim_verbose >= LOG_SOLVER)
+    {
+      sim_result->emit();
+    }
+  //saveall();
+  //update_DAEsystem();
+
+  storeExtrapolationData();
+  storeExtrapolationData();
+
+  if (sim_verbose >= LOG_SOLVER)
+    {
+      sim_result->emit();
+    }
+
+  //first try with the given method as default simplex and
+  //then try with the other one
+  int retVal = 0;
+  retVal = initialize(init_method);
+  if (retVal != 0)
+    {
+      if (init_method == std::string("simplex"))
+        {
+          init_method = std::string("newuoa");
+          retVal = initialize(init_method);
+        }
+      else if (init_method == std::string("newuoa"))
+        {
+          init_method = std::string("simplex");
+          retVal = initialize(init_method);
+        }
+      if (retVal != 0)
+        {
+          printf("Initialization of the current initial set of equations and initial guesses fails!\n");
+          printf("Try with better Initial guesses for the states.\n");
+          retVal = 1;
+        }
+    }
+
+  saveall();
+  storeExtrapolationData();
+  storeExtrapolationData();
+
+  if (sim_verbose >= LOG_SOLVER)
+    {
+      sim_result->emit();
+    }
+
+  update_DAEsystem();
+  SaveZeroCrossings();
+  saveall();
+  if (sim_verbose >= LOG_SOLVER)
+    {
+      sim_result->emit();
+    }
+  storeExtrapolationData();
+  storeExtrapolationData();
+  globalData->init = 0;
+  return retVal;
+}
+

c_runtime/simulation_init.h

@@ -39,7 +39,7 @@
 #ifndef _SIMULATION_INIT_H
 #define _SIMULATION_INIT_H
 
-int initialize(const std::string*method);
+int main_initialize(const std::string*method);
 
 #ifndef NEWUOA
 #define NEWUOA newuoa_

c_runtime/simulation_input.cpp

@@ -132,66 +132,66 @@ void read_commented_value(ifstream &f, signed char *str);
   }
   for(int i = 0; i < simData->nStates; i++) { // Read x initial values
     read_commented_value(file,&simData->states[i]);
-    if (sim_verbose >= LOG_SOLVER) {
+    if (sim_verbose >= LOG_INIT) {
       cout << "read " << simData->statesNames[i].name << " = " << simData->states[i] << " from init file." << endl;
     }
   }
   for(int i = 0; i < simData->nStates; i++) { // Read der(x) initial values
     read_commented_value(file,&simData->statesDerivatives[i]);
-    if (sim_verbose >= LOG_SOLVER) {
+    if (sim_verbose >= LOG_INIT) {
       cout << "read " << simData->stateDerivativesNames[i].name << " = " << simData->statesDerivatives[i] << " from init file." << endl;
     }
   }
   for(int i = 0; i < simData->nAlgebraic; i++) { // Read y initial values
     read_commented_value(file,&simData->algebraics[i]);
-    if (sim_verbose >= LOG_SOLVER) {
+    if (sim_verbose >= LOG_INIT) {
       cout << "read " << simData->algebraicsNames[i].name << " = " << simData->algebraics[i] << " from init file." << endl;
     }
   }
   for(int i = 0; i < simData->nParameters; i++) { // Read parameter values
     read_commented_value(file,&simData->parameters[i]);
-    if (sim_verbose >= LOG_SOLVER) {
+    if (sim_verbose >= LOG_INIT) {
       cout << "read " << simData->parametersNames[i].name << " = " << simData->parameters[i] << " from init file." << endl;
     }
   }
 
   for(int i = 0; i < simData->intVariables.nParameters; i++) { // Read parameter values
     read_commented_value(file,&simData->intVariables.parameters[i]);
-    if (sim_verbose >= LOG_SOLVER) {
+    if (sim_verbose >= LOG_INIT) {
       cout << "read " << simData->int_param_names[i].name << " = " << simData->intVariables.parameters[i] << " from init file." << endl;
     }
   }
 
   for(int i = 0; i < simData->intVariables.nAlgebraic; i++) { // Read parameter values
     read_commented_value(file,&simData->intVariables.algebraics[i]);
-    if (sim_verbose >= LOG_SOLVER) {
+    if (sim_verbose >= LOG_INIT) {
       cout << "read " << simData->int_alg_names[i].name << " = " << simData->intVariables.algebraics[i] << " from init file." << endl;
     }
   }
 
   for(int i = 0; i < simData->boolVariables.nParameters; i++) { // Read parameter values
     read_commented_value(file,&simData->boolVariables.parameters[i]);
-    if (sim_verbose >= LOG_SOLVER) {
+    if (sim_verbose >= LOG_INIT) {
       cout << "read " << simData->bool_param_names[i].name << " = " << (bool)simData->boolVariables.parameters[i] << " from init file." << endl;
     }
   }
 
   for(int i = 0; i < simData->boolVariables.nAlgebraic; i++) { // Read parameter values
     read_commented_value(file,&simData->boolVariables.algebraics[i]);
-    if (sim_verbose >= LOG_SOLVER) {
+    if (sim_verbose >= LOG_INIT) {
       cout << "read " << simData->bool_alg_names[i].name << " = " << (bool)simData->boolVariables.algebraics[i] << " from init file." << endl;
     }
   }
 
   for(int i=0; i < simData->stringVariables.nParameters; i++) { // Read string parameter values
     read_commented_value(file,&(simData->stringVariables.parameters[i]));
-    if (sim_verbose >= LOG_SOLVER) {
+    if (sim_verbose >= LOG_INIT) {
       cout << "read " << simData->string_param_names[i].name << " = \"" << simData->stringVariables.parameters[i] << "\" from init file." << endl;
     }
   }
   for(int i=0; i < simData->stringVariables.nAlgebraic; i++) { // Read string algebraic values
     read_commented_value(file,&simData->stringVariables.algebraics[i]);
-    if (sim_verbose >= LOG_SOLVER) {
+    if (sim_verbose >= LOG_INIT) {
       cout << "read " << simData->string_alg_names[i].name << " from init file." << endl;
     }
   }

c_runtime/simulation_runtime.cpp

@@ -100,6 +100,7 @@ simulation_result *sim_result;
 /* Flags for controlling logging to stdout */
 const int LOG_STATS = 1;
 const int LOG_INIT = 2;
+const int LOG_RES_INIT = 3;
 const int LOG_SOLVER = 4;
 const int LOG_NONLIN_SYS = 8;
 const int LOG_EVENTS = 16;
@@ -232,10 +233,14 @@ verboseLevel(int argc, char**argv)
     {
       res |= LOG_STATS;
     }
-if (flags->find("LOG_INIT", 0) != string::npos)
+  if (flags->find("LOG_INIT", 0) != string::npos)
     {
       res |= LOG_INIT;
     }
+  if (flags->find("LOG_RES_INIT", 0) != string::npos)
+    {
+      res |= LOG_RES_INIT;
+    }
   if (flags->find("LOG_SOLVER", 0) != string::npos)
     {
       res |= LOG_SOLVER;

c_runtime/simulation_runtime.h

@@ -103,6 +103,7 @@ extern omc_fileInfo TermInfo; /* message for termination. */
 /* Flags for controlling logging to stdout */
 extern const int LOG_STATS;
 extern const int LOG_INIT;
+extern const int LOG_RES_INIT;
 extern const int LOG_SOLVER;
 extern const int LOG_EVENTS;
 extern const int LOG_NONLIN_SYS;
@@ -223,6 +224,7 @@ typedef struct sim_DATA {
   double* statesBackup;
 
   char* initFixed; /* Fixed attribute for all variables and parameters */
+  char* var_attr; /* Type attribute for all variables and parameters */
   int init; /* =1 during initialization, 0 otherwise. */
   void** extObjs; /* External objects */
   /* nStatesDerivatives == states */