Enable global homotopy as fallback option by default

Introduce runtime flag homotopyOnFirstTry (default: false).
If that flag is not activated first try to solve the initialization
problem without homotopy.

Author: ptaeuber

Compiler/SimCode/SimCodeUtil.mo

@@ -261,13 +261,22 @@ algorithm
     end if;
 
     // initialization stuff
-    (initialEquations, removedInitialEquations, uniqueEqIndex, tempvars) := createInitialEquations(inInitDAE, inRemovedInitialEquationLst, uniqueEqIndex, {});
+    // ********************
+
+    // generate equations for initDAE
+    (initialEquations, uniqueEqIndex, tempvars) := createInitialEquations(inInitDAE, uniqueEqIndex, {});
+
+    // generate equations for initDAE_lambda0
     if isSome(inInitDAE_lambda0) then
       SOME(initDAE_lambda0) := inInitDAE_lambda0;
       (initialEquations_lambda0, uniqueEqIndex, tempvars) := createInitialEquations_lambda0(initDAE_lambda0, uniqueEqIndex, tempvars);
     else
       initialEquations_lambda0 := {};
     end if;
+
+    // generate equations for removed initial equations
+    (removedInitialEquations, uniqueEqIndex, tempvars) := createNonlinearResidualEquations(inRemovedInitialEquationLst, uniqueEqIndex, tempvars);
+
     execStat("simCode: created initialization part");
 
     shared as BackendDAE.SHARED(globalKnownVars=globalKnownVars,
@@ -6358,18 +6367,16 @@ end createSingleAlgorithmCode;
 
 protected function createInitialEquations "author: lochel"
   input BackendDAE.BackendDAE inInitDAE;
-  input List<BackendDAE.Equation> inRemovedEqnLst;
   input Integer iuniqueEqIndex;
   input list<SimCodeVar.SimVar> itempvars;
   output list<SimCode.SimEqSystem> outInitialEqns = {};
-  output list<SimCode.SimEqSystem> outRemovedInitialEqns = {};
   output Integer ouniqueEqIndex = iuniqueEqIndex;
   output list<SimCodeVar.SimVar> otempvars = itempvars;
 protected
   BackendDAE.EquationArray  removedEqs;
   list<SimCodeVar.SimVar> tempvars;
   Integer uniqueEqIndex;
-  list<SimCode.SimEqSystem> allEquations, knownVarEquations, solvedEquations, removedEquations, aliasEquations, removedInitialEquations;
+  list<SimCode.SimEqSystem> allEquations, knownVarEquations, solvedEquations, removedEquations, aliasEquations;
   BackendDAE.EqSystems systs;
   BackendDAE.Shared shared;
   BackendDAE.Variables globalKnownVars, aliasVars;
@@ -6391,12 +6398,8 @@ algorithm
   allEquations := List.append_reverse(solvedEquations, allEquations);
   allEquations := listAppend(knownVarEquations, allEquations);
 
-  // generate equations from removed initial equations
-  (removedInitialEquations, uniqueEqIndex, tempvars) := createNonlinearResidualEquations(inRemovedEqnLst, uniqueEqIndex, tempvars);
-
   // output
   outInitialEqns := allEquations;
-  outRemovedInitialEqns := removedInitialEquations;
   ouniqueEqIndex := uniqueEqIndex;
   otempvars := tempvars;
 end createInitialEquations;

Compiler/Util/Flags.mo

@@ -1373,10 +1373,10 @@ constant ConfigFlag ZEROMQ_FILE_SUFFIX = CONFIG_FLAG(115, "zeroMQFileSuffix",
   SOME("z"), EXTERNAL(), STRING_FLAG(""), NONE(),
   Util.gettext("Sets the file suffix for zeroMQ port file if --interactive=zmq is used."));
 constant ConfigFlag HOMOTOPY_APPROACH = CONFIG_FLAG(116, "homotopyApproach",
-  NONE(), EXTERNAL(), STRING_FLAG("local"),
+  NONE(), EXTERNAL(), STRING_FLAG("global"),
   SOME(STRING_DESC_OPTION({
-    ("local", Util.gettext("Default, local homotopy approach. The homotopy parameter only effects the local strongly connected component.")),
-    ("global", Util.gettext("Global homotopy approach. The homotopy parameter effects the entire initialization system."))
+    ("local", Util.gettext("Local homotopy approach. The homotopy parameter only effects the local strongly connected component.")),
+    ("global", Util.gettext("Default, global homotopy approach. The homotopy parameter effects the entire initialization system."))
     })),
     Util.gettext("Sets the homotopy approach."));
 

SimulationRuntime/c/simulation/solver/initialization/initialization.c

@@ -37,6 +37,7 @@
 #include "util/omc_error.h"
 #include "openmodelica.h"
 #include "openmodelica_func.h"
+#include "simulation/options.h"
 #include "simulation/solver/model_help.h"
 #if !defined(OMC_MINIMAL_RUNTIME)
 #include "util/read_matlab4.h"
@@ -73,7 +74,7 @@
 #include <math.h>
 #include <string.h>
 
-int init_lambda_steps = 1;
+int init_lambda_steps = 4;
 
 /*! \fn void dumpInitializationStatus(DATA *data)
  *
@@ -193,6 +194,39 @@ static int symbolic_initialization(DATA *data, threadData_t *threadData)
   storePreValues(data);
   overwriteOldSimulationData(data);
 
+  /* If there is no homotopy in the model or local homotopy is activated
+     or homotopy is disabled by runtime flag '-ils=<lambda_steps>',
+     solve WITHOUT HOMOTOPY. */
+  if (data->callback->useHomotopy == 0 || init_lambda_steps < 2){
+    data->simulationInfo->lambda = 1.0;
+    data->callback->functionInitialEquations(data, threadData);
+
+  /* If there is homotopy in the model and global homotopy is activated
+     and homotopy on first try is deactivated,
+     TRY TO SOLVE WITHOUT HOMOTOPY FIRST. */
+  } else if (!omc_flag[FLAG_HOMOTOPY_ON_FIRST_TRY]) {
+    /* try */
+#ifndef OMC_EMCC
+  MMC_TRY_INTERNAL(simulationJumpBuffer)
+#endif
+
+    data->simulationInfo->lambda = 1.0;
+    infoStreamPrint(LOG_INIT, 0, "Try to solve the initialization problem without homotopy first.");
+    data->callback->functionInitialEquations(data, threadData);
+    init_lambda_steps = 0;
+    infoStreamPrint(LOG_INIT, 0, "Initialization finished without homotopy.");
+
+    /* catch */
+#ifndef OMC_EMCC
+  MMC_CATCH_INTERNAL(simulationJumpBuffer)
+#endif
+    if(init_lambda_steps > 0)
+      warningStreamPrint(LOG_ASSERT, 0, "Failed to solve the initial system without homotopy. If homotopy is available the homotopy method is now used.");
+  }
+
+  /* If there is homotopy in the model and global homotopy is activated
+     and solving without homotopy failed or is not wanted,
+     use GLOBAL HOMOTOPY METHOD. */
   if (data->callback->useHomotopy == 1 && init_lambda_steps > 1)
   {
     long step;
@@ -203,15 +237,15 @@ static int symbolic_initialization(DATA *data, threadData_t *threadData)
     {
       sprintf(buffer, "%s_global_homotopy.csv", mData->modelFilePrefix);
       pFile = fopen(buffer, "wt");
-      fprintf(pFile, "%s", "lambda");
+      fprintf(pFile, "\"sep=,\"\n%s", "lambda");
       for(i=0; i<mData->nVariablesReal; ++i)
         fprintf(pFile, ",%s", mData->realVarsData[i].info.name);
       fprintf(pFile, "\n");
     }
 #endif
 
     infoStreamPrint(LOG_INIT, 1, "homotopy process\n---------------------------");
-    for(step=0; step<init_lambda_steps-1; ++step)
+    for(step=0; step<init_lambda_steps; ++step)
     {
       data->simulationInfo->lambda = ((double)step)/(init_lambda_steps-1);
       infoStreamPrint(LOG_INIT, 0, "homotopy parameter lambda = %g", data->simulationInfo->lambda);
@@ -243,23 +277,13 @@ static int symbolic_initialization(DATA *data, threadData_t *threadData)
         break;
     }
     messageClose(LOG_INIT);
-    infoStreamPrint(LOG_INIT, 0, "homotopy parameter lambda = 1");
   }
 
-  data->simulationInfo->lambda = 1.0;
-  data->callback->functionInitialEquations(data, threadData);
   storeRelations(data);
 
 #if !defined(OMC_NO_FILESYSTEM)
   if(data->callback->useHomotopy == 1 && init_lambda_steps > 1 && ACTIVE_STREAM(LOG_INIT))
-  {
-    infoStreamPrint(LOG_INIT, 0, "homotopy parameter lambda = %g done\n---------------------------", data->simulationInfo->lambda);
-    fprintf(pFile, "%.16g", data->simulationInfo->lambda);
-    for(i=0; i<mData->nVariablesReal; ++i)
-      fprintf(pFile, ",%.16g", data->localData[0]->realVars[i]);
-    fprintf(pFile, "\n");
     fclose(pFile);
-  }
 #endif
 
   /* check for over-determined systems */

SimulationRuntime/c/util/simulation_options.c

@@ -46,6 +46,7 @@ const char *FLAG_NAME[FLAG_MAX+1] = {
   /* FLAG_EMIT_PROTECTED */        "emit_protected",
   /* FLAG_F */                     "f",
   /* FLAG_HELP */                  "help",
+  /* FLAG_HOMOTOPY_ON_FIRST_TRY */ "homotopyOnFirstTry",
   /* FLAG_IDA_MAXERRORTESTFAIL */  "idaMaxErrorTestFails",
   /* FLAG_IDA_MAXNONLINITERS */    "idaMaxNonLinIters",
   /* FLAG_IDA_MAXCONVFAILS */      "idaMaxConvFails",
@@ -134,6 +135,7 @@ const char *FLAG_DESC[FLAG_MAX+1] = {
   /* FLAG_EMIT_PROTECTED */        "emits protected variables to the result-file",
   /* FLAG_F */                     "value specifies a new setup XML file to the generated simulation code",
   /* FLAG_HELP */                  "get detailed information that specifies the command-line flag",
+  /* FLAG_HOMOTOPY_ON_FIRST_TRY */ "Directly use the homotopy method to solve the initialization problem.",
   /* FLAG_IDA_MAXERRORTESTFAIL */  "value specifies the maximum number of error test failures in attempting one step. The default value is 7.",
   /* FLAG_IDA_MAXNONLINITERS */    "value specifies the maximum number of nonlinear solver iterations at one step. The default value is 3.",
   /* FLAG_IDA_MAXCONVFAILS */      "value specifies the maximum number of nonlinear solver convergence failures at one step. The default value is 10.",
@@ -144,7 +146,7 @@ const char *FLAG_DESC[FLAG_MAX+1] = {
   /* FLAG_IIF */                   "value specifies an external file for the initialization of the model",
   /* FLAG_IIM */                   "value specifies the initialization method",
   /* FLAG_IIT */                   "[double] value specifies a time for the initialization of the model",
-  /* FLAG_ILS */                   "[int] default: 1",
+  /* FLAG_ILS */                   "[int] default: 4",
   /* FLAG_IMPRK_ORDER */           "[int (default 5)] value specifies the integration order of the implicit Runge-Kutta method. Valid values: 1-6",
   /* FLAG_IMPRK_LS */              "selects the linear solver of the integration methods: impeuler, trapezoid and imprungekuta",
   /* FLAG_INITIAL_STEP_SIZE */     "value specifies an initial step size for supported solver",
@@ -242,6 +244,9 @@ const char *FLAG_DETAILED_DESC[FLAG_MAX+1] = {
   /* FLAG_HELP */
   "  Get detailed information that specifies the command-line flag\n"
   "  For example, -help=f prints detailed information for command-line flag f.",
+  /* FLAG_HOMOTOPY_ON_FIRST_TRY */
+  "  If the model contains the homotopy operator, directly use the homotopy method to solve the initialization problem.\n"
+  "  If disabled, first try to solve without homotopy and only use homotopy as fallback option.",
   /* FLAG_IDA_MAXERRORTESTFAIL */
   "  value specifies the maximum number of error test failures in attempting one step. The default value is 7.",
   /* FLAG_IDA_MAXNONLINITERS */
@@ -269,7 +274,7 @@ const char *FLAG_DETAILED_DESC[FLAG_MAX+1] = {
   "  Value [Real] specifies a time for the initialization of the model.",
   /* FLAG_ILS */
   "  Value specifies the number of steps for homotopy method (required: -iim=symbolic).\n"
-  "  The value is an Integer with default value 1.",
+  "  The value is an Integer with default value 4.",
   /* FLAG_IMPRK_ORDER */
   "  Value specifies the integration order of the implicit Runge-Kutta method. Valid values: 1 to 6. Default order is 5.",
   /* FLAG_IMPRK_LS */
@@ -452,6 +457,7 @@ const int FLAG_TYPE[FLAG_MAX] = {
   /* FLAG_EMIT_PROTECTED */        FLAG_TYPE_FLAG,
   /* FLAG_F */                     FLAG_TYPE_OPTION,
   /* FLAG_HELP */                  FLAG_TYPE_OPTION,
+  /* FLAG_HOMOTOPY_ON_FIRST_TRY */ FLAG_TYPE_FLAG,
   /* FLAG_IDA_MAXERRORTESTFAIL */  FLAG_TYPE_OPTION,
   /* FLAG_IDA_MAXNONLINITERS */    FLAG_TYPE_OPTION,
   /* FLAG_IDA_MAXCONVFAILS */      FLAG_TYPE_OPTION,

SimulationRuntime/c/util/simulation_options.h

@@ -54,6 +54,7 @@ enum _FLAG
   FLAG_EMIT_PROTECTED,
   FLAG_F,
   FLAG_HELP,
+  FLAG_HOMOTOPY_ON_FIRST_TRY,
   FLAG_IDA_MAXERRORTESTFAIL,
   FLAG_IDA_MAXNONLINITERS,
   FLAG_IDA_MAXCONVFAILS,