Initialization with homotopy and --daeMode does not generate initial-lambda0 equations

[Francesco Casella]

When the homotopy operator is used, a separate set of intial equations, labelled as initial-lambda0 in the equations browser of the transformational debugger, is generated to solve the simplified system corresponding to lambda = 0, where lambda is the homotopy parameter. This is done in order to exploit the structural simplifications induced by the simplified expressions in the model, which help achieve convergence at the initialization of the homotopy process.

Unfortunately this step is skipped when --daeMode is activated. This is problematic for all the power grids models that need to be solved with that strategy. The attached model demonstrates the issue.

I may be unaware of some technical implementation details, but as far as I understand, the initialization phase is exactly the same whether the dynamic equations are solved with odeMode or daeMode. I guess when daeMode was implemented Willi didn't care about this aspect, but porting it from odeMode should be easy.

Could you please do that?

Thank you!

[Francesco Casella]

See the ​PowerGrids CI report, in particular ​IEEE14busLoad2Variation

[Karim Adbdelhak]

I added the lambda0 equation system in ​PR728 to DAEMode.

If you have any benchmarks to test this please do so! :)

[Francesco Casella]

Thank you @Karim, I will test the nighly build tomorrow morning and report.

[Francesco Casella]

@Karim, your PR broke several models of PowerGrids that use --daeMode, see ​report. I hope it's just due to some minor oversight.

I would suggest that you revert your PR for the time being so we restore the previous behaviour on the nightlies, and then test it locally using those tests from the library that you can find on ​https://github.com/powergrids.

I cannot download the nightly build right now, so I cannot get any further in the analysis. I'll try to do that later on.

[Francesco Casella]

Retargeted to 1.17.0 after 1.16.0 release

[Francesco Casella]

Related to #6338.

I'm not sure what the current status is. Karim first pushed ​PR 728, but that broke ​several test cases in PowerGrids. Then I asked Karim to revert the push, which he did in ​PR 732. I'm not sure if there has been any progress.

I would suggest to restart from ​PR 728 and figure out what's wrong there. I can help if necessary, this is quite critical for my talk at the Workshop.

[Andreas Heuermann]

I'm working with this test script at the moment to see why PR 728 was failing for models from PowerGrids:

cd("Temp"); getErrorString();
loadFile("/home/andreas/workspace/OpenModelica/build/lib/omlibrary/PowerGrids 1.0.0/package.mo"); getErrorString();
setCommandLineOptions("-d=newInst,dumpSimCode --daeMode"); getErrorString();
simulate(PowerGrids.Examples.Tutorial.GridOperation.Controlled.ControlledGridParametersChange); getErrorString();

Debugging shows that we seem to have some kind of problem with de-initializing nonlinear systems in ​freeNonlinearSystems() while calling

free(nonlinsys[i].nlsxOld);

for i=0.
I'm getting
free(): invalid next size (fast).

The reason: We are allocating nlsxOld with space for 4 doubles:

nonlinsys[0].nlsxOld = (double*) malloc(4*sizeof(double));

but generate code that access it at other places:

  /* get old value */
  data->simulationInfo->nonlinearSystemData[0].nlsxOld[0] = data->simulationInfo->realParameter[81] /* GEN.iqPuStart PARAM */;
  data->simulationInfo->nonlinearSystemData[0].nlsxOld[1] = data->simulationInfo->realParameter[83] /* GEN.lambdadPuStart PARAM */;
  data->simulationInfo->nonlinearSystemData[0].nlsxOld[2] = data->simulationInfo->realParameter[80] /* GEN.ifPuStart PARAM */;
  data->simulationInfo->nonlinearSystemData[0].nlsxOld[3] = data->simulationInfo->realParameter[112] /* GEN.ufPuStart PARAM */;
  data->simulationInfo->nonlinearSystemData[0].nlsxOld[4] = data->simulationInfo->realParameter[111] /* GEN.ufPuInStart PARAM */;
  data->simulationInfo->nonlinearSystemData[0].nlsxOld[5] = data->localData[0]->realVars[34] /* AVR.firstOrderLim.Gk.y variable */;
  data->simulationInfo->nonlinearSystemData[0].nlsxOld[6] = data->localData[0]->realVars[42] /* AVR.variableLimiter.y variable */;
  data->simulationInfo->nonlinearSystemData[0].nlsxOld[7] = data->simulationInfo->realParameter[35] /* AVR.leadLag.yStart PARAM */;
  data->simulationInfo->nonlinearSystemData[0].nlsxOld[8] = data->localData[0]->realVars[40] /* AVR.leadLag.feedback.y variable */;
  data->simulationInfo->nonlinearSystemData[0].nlsxOld[9] = data->localData[0]->realVars[39] /* AVR.inputLimiter.y variable */;
  data->simulationInfo->nonlinearSystemData[0].nlsxOld[10] = data->localData[0]->realVars[38] /* AVR.inputLimiter.u variable */;
  data->simulationInfo->nonlinearSystemData[0].nlsxOld[11] = DIVISION_SIM(data->simulationInfo->realParameter[98] /* GEN.port.VStart PARAM */,12124.35565298214,"12124.35565298214",equationIndexes);
  data->simulationInfo->nonlinearSystemData[0].nlsxOld[12] = data->simulationInfo->realParameter[96] /* GEN.port.UStart PARAM */;
  data->simulationInfo->nonlinearSystemData[0].nlsxOld[13] = data->localData[0]->realVars[70] /* GEN.port.u.im variable */;
  data->simulationInfo->nonlinearSystemData[0].nlsxOld[14] = data->localData[0]->realVars[249] /* TGEN.vA.im variable */;
  data->simulationInfo->nonlinearSystemData[0].nlsxOld[15] = data->localData[0]->realVars[251] /* TGEN.vAt.im variable */;
  data->simulationInfo->nonlinearSystemData[0].nlsxOld[16] = data->localData[0]->realVars[253] /* TGEN.vz.im variable */;
  data->simulationInfo->nonlinearSystemData[0].nlsxOld[17] = data->localData[0]->realVars[207] /* TGEN.iAt.im variable */;
  data->simulationInfo->nonlinearSystemData[0].nlsxOld[18] = data->localData[0]->realVars[208] /* TGEN.iAt.re variable */;
  data->simulationInfo->nonlinearSystemData[0].nlsxOld[19] = data->localData[0]->realVars[204] /* TGEN.iA.re variable */;
  data->simulationInfo->nonlinearSystemData[0].nlsxOld[20] = data->localData[0]->realVars[69] /* GEN.port.iGenPu.re variable */;
  data->simulationInfo->nonlinearSystemData[0].nlsxOld[21] = data->simulationInfo->realParameter[109] /* GEN.thetaStart PARAM */;
  data->simulationInfo->nonlinearSystemData[0].nlsxOld[22] = data->simulationInfo->realParameter[110] /* GEN.udPuStart PARAM */;
  data->simulationInfo->nonlinearSystemData[0].nlsxOld[23] = data->simulationInfo->realParameter[84] /* GEN.lambdaqPuStart PARAM */;
  data->simulationInfo->nonlinearSystemData[0].nlsxOld[24] = data->simulationInfo->realParameter[79] /* GEN.idPuStart PARAM */;
  data->simulationInfo->nonlinearSystemData[0].nlsxOld[25] = data->localData[0]->realVars[68] /* GEN.port.iGenPu.im variable */;
  data->simulationInfo->nonlinearSystemData[0].nlsxOld[26] = data->localData[0]->realVars[203] /* TGEN.iA.im variable */;
  data->simulationInfo->nonlinearSystemData[0].nlsxOld[27] = data->simulationInfo->realParameter[113] /* GEN.uqPuStart PARAM */;
  data->simulationInfo->nonlinearSystemData[0].nlsxOld[28] = data->localData[0]->realVars[254] /* TGEN.vz.re variable */;
  data->simulationInfo->nonlinearSystemData[0].nlsxOld[29] = data->localData[0]->realVars[252] /* TGEN.vAt.re variable */;
  data->simulationInfo->nonlinearSystemData[0].nlsxOld[30] = data->localData[0]->realVars[250] /* TGEN.vA.re variable */;
  data->simulationInfo->nonlinearSystemData[0].nlsxOld[31] = data->localData[0]->realVars[71] /* GEN.port.u.re variable */;
  data->simulationInfo->nonlinearSystemData[0].nlsxOld[32] = data->localData[0]->realVars[73] /* GEN.port.vPu.re variable */;
  data->simulationInfo->nonlinearSystemData[0].nlsxOld[33] = data->localData[0]->realVars[72] /* GEN.port.vPu.im variable */;
  data->simulationInfo->nonlinearSystemData[0].nlsxOld[34] = data->simulationInfo->realParameter[32] /* AVR.leadLag.firstOrder.y_start PARAM */;

In that location it should probably be nonlinearSystemData[5].

[Andreas Heuermann]

@kabdelhak and myself found the solution.

See ​PR#7100.

[Francesco Casella]

I love you, great team! :)

I just started a Windows nightly, will check later today with my test. Did you try it already?

[Andreas Heuermann]

The example above works fine, but I didn't look at any results and I didn't test other examples.

[Francesco Casella]

Type1_N_4_M_4 now initializes successfully in a reasonable amount of time, given that it involves three homotopy steps, about 100 ms. I will now continue with larger tests, and report during the Workshop.

Thank you guys, I appreciate this a lot!