Skip to content

Commit

Permalink
[C/DO]: fixed values for discrete variables stored in the section for…
Browse files Browse the repository at this point in the history
… real variables (#8036)

* fix: issue4979

* refacator code avoid duplicated code
* algorithmen, handling only state and inputs
  * copy of the real discret value was missing

* test: update output

* TFC5 was broken.
  * updated the initial guess.
  • Loading branch information
vruge committed Oct 22, 2021
1 parent 70dd572 commit 998ecd7
Show file tree
Hide file tree
Showing 12 changed files with 280 additions and 143 deletions.
Expand Up @@ -110,7 +110,7 @@ void debugeJac(OptData * optData, Number* vopt){
const int npv = np*nv;
const int nt = optData->dim.nt;
const int NRes = optData->dim.NRes;
const int nReal = optData->data->modelData->nVariablesReal;
const int nReal = optData->dim.nReal;
const int NV = optData->dim.NV;
Number vopt_shift[NV];
long double h[nv][nsi][np];
Expand Down
Expand Up @@ -142,9 +142,8 @@ inline void allocate_der_struct(OptDataStructure *s, OptDataDim * dim, DATA* dat
}

if(s->mayer){
const int nReal = dim->nReal;
dim->index_mayer = -1;
for(i = 0; i < nReal; ++i){
for(i = 0; i < dim->nReal; ++i){
if(&data->localData[0]->realVars[i] == s->pmayer){
dim->index_mayer = i;
break;
Expand All @@ -153,9 +152,8 @@ inline void allocate_der_struct(OptDataStructure *s, OptDataDim * dim, DATA* dat
}

if(s->lagrange){
const int nReal = dim->nReal;
dim->index_lagrange = -1;
for(i = 0; i < nReal; ++i){
for(i = 0; i < dim->nReal; ++i){
if(&data->localData[0]->realVars[i] == s->plagrange){
dim->index_lagrange = i;
break;
Expand Down
Expand Up @@ -136,9 +136,6 @@ static short initial_guess_ipopt_sim(OptData *optData, SOLVER_INFO* solverInfo,

if((double)data->simulationInfo->startTime < optData->time.t0){
double t = data->simulationInfo->startTime;
const int nBoolean = data->modelData->nVariablesBoolean;
const int nInteger = data->modelData->nVariablesInteger;
const int nRelations = data->modelData->nRelations;

FILE * pFile = optData->pFile;
fprintf(pFile, "%lf ",(double)t);
Expand All @@ -162,15 +159,7 @@ static short initial_guess_ipopt_sim(OptData *optData, SOLVER_INFO* solverInfo,
}
fprintf(pFile, "%s", "\n");
}
memcpy(optData->v0, data->localData[0]->realVars, nReal*sizeof(modelica_real));
memcpy(optData->i0, data->localData[0]->integerVars, nInteger*sizeof(modelica_integer));
memcpy(optData->b0, data->localData[0]->booleanVars, nBoolean*sizeof(modelica_boolean));
memcpy(optData->i0Pre, data->simulationInfo->integerVarsPre, nInteger*sizeof(modelica_integer));
memcpy(optData->b0Pre, data->simulationInfo->booleanVarsPre, nBoolean*sizeof(modelica_boolean));
memcpy(optData->v0Pre, data->simulationInfo->realVarsPre, nReal*sizeof(modelica_real));
memcpy(optData->rePre, data->simulationInfo->relationsPre, nRelations*sizeof(modelica_boolean));
memcpy(optData->re, data->simulationInfo->relations, nRelations*sizeof(modelica_boolean));
memcpy(optData->storeR, data->simulationInfo->storedRelations, nRelations*sizeof(modelica_boolean));
copy_initial_values(optData, data);

if(1){
printf("\n--------------------------------------------------------");
Expand Down
Expand Up @@ -54,8 +54,7 @@ static inline void pickUpStates(OptData* optdata);
static inline void updateDOSystem(OptData * optData, DATA * data, threadData_t *threadData,
const int i, const int j, const int index, const int m);

static inline void setLocalVars(OptData * optData, DATA * data, const double * const vopt,
const int i, const int j, const int shift);
void setLocalVars(OptData * optData, DATA * data, const double * const vopt, const int i, const int j, const int shift);

static inline int getNsi(char*, const int, modelica_boolean*);
static inline void overwriteTimeGridFile(OptDataTime * time, char* filename, long double c[], const int np, const int nsi);
Expand Down Expand Up @@ -606,7 +605,6 @@ void res2file(OptData *optData, SOLVER_INFO* solverInfo, double *vopt){
SIMULATION_DATA *sData = (SIMULATION_DATA*)data->localData[0];

FILE * pFile = optData->pFile;
double * v0 = optData->v0;
double *vnom = optData->bounds.vnom;
long double **t = optData->time.t;
long double t0 = optData->time.t0;
Expand Down Expand Up @@ -638,15 +636,7 @@ void res2file(OptData *optData, SOLVER_INFO* solverInfo, double *vopt){
}
fprintf(pFile, "%s", "\n");
/******************/
memcpy(sData->realVars, v0, nReal*sizeof(modelica_real));
memcpy(data->localData[0]->integerVars, optData->i0, nInteger*sizeof(modelica_integer));
memcpy(data->localData[0]->booleanVars, optData->b0, nBoolean*sizeof(modelica_boolean));
memcpy(data->simulationInfo->integerVarsPre, optData->i0Pre, nInteger*sizeof(modelica_integer));
memcpy(data->simulationInfo->booleanVarsPre, optData->b0Pre, nBoolean*sizeof(modelica_boolean));
memcpy(data->simulationInfo->realVarsPre, optData->v0Pre, nReal*sizeof(modelica_real));
memcpy(data->simulationInfo->relationsPre, optData->rePre, nRelations*sizeof(modelica_boolean));
memcpy(data->simulationInfo->relations, optData->re, nRelations*sizeof(modelica_boolean));
memcpy(data->simulationInfo->storedRelations, optData->storeR, nRelations*sizeof(modelica_boolean));
copy_initial_values(optData, data);
/******************/
solverInfo->currentTime = (double)t0;
sData->timeValue = solverInfo->currentTime;
Expand Down Expand Up @@ -679,6 +669,25 @@ void res2file(OptData *optData, SOLVER_INFO* solverInfo, double *vopt){
fclose(pFile);
}


void copy_initial_values(OptData * optData, DATA* data){
const int nBoolean = optData->data->modelData->nVariablesBoolean ;
const int nInteger = optData->data->modelData->nVariablesInteger;
const int nReal = optData->dim.nReal;
const int nRelations = optData->data->modelData->nRelations;

memcpy(data->localData[0]->realVars, optData->v0, nReal*sizeof(modelica_real));
memcpy(data->localData[0]->integerVars, optData->i0, nInteger*sizeof(modelica_integer));
memcpy(data->localData[0]->booleanVars, optData->b0, nBoolean*sizeof(modelica_boolean));
memcpy(data->simulationInfo->integerVarsPre, optData->i0Pre, nInteger*sizeof(modelica_integer));
memcpy(data->simulationInfo->booleanVarsPre, optData->b0Pre, nBoolean*sizeof(modelica_boolean));
memcpy(data->simulationInfo->realVarsPre, optData->v0Pre, nReal*sizeof(modelica_real));
memcpy(data->simulationInfo->relationsPre, optData->rePre, nRelations*sizeof(modelica_boolean));
memcpy(data->simulationInfo->relations, optData->re, nRelations*sizeof(modelica_boolean));
memcpy(data->simulationInfo->storedRelations, optData->storeR, nRelations*sizeof(modelica_boolean));

}

/*!
* transfer optimizer data to model data
* author: Vitalij Ruge
Expand All @@ -688,12 +697,6 @@ void optData2ModelData(OptData *optData, double *vopt, const int index){
const int nsi = optData->dim.nsi;
const int np = optData->dim.np;

const int nBoolean = optData->data->modelData->nVariablesBoolean;
const int nInteger = optData->data->modelData->nVariablesInteger;
const int nReal = optData->dim.nReal;
const int nRelations = optData->data->modelData->nRelations;


modelica_real * realVars[3];
modelica_real * tmpVars[2] = {NULL, NULL};

Expand All @@ -709,16 +712,7 @@ void optData2ModelData(OptData *optData, double *vopt, const int index){
if(optData->s.matrix[l])
tmpVars[l] = data->simulationInfo->analyticJacobians[indexBC[l]].tmpVars;
}

memcpy(data->localData[0]->integerVars, optData->i0, nInteger*sizeof(modelica_integer));
memcpy(data->localData[0]->booleanVars, optData->b0, nBoolean*sizeof(modelica_boolean));
memcpy(data->simulationInfo->integerVarsPre, optData->i0Pre, nInteger*sizeof(modelica_integer));
memcpy(data->simulationInfo->booleanVarsPre, optData->b0Pre, nBoolean*sizeof(modelica_boolean));
memcpy(data->simulationInfo->realVarsPre, optData->v0Pre, nReal*sizeof(modelica_real));
memcpy(data->simulationInfo->relationsPre, optData->rePre, nRelations*sizeof(modelica_boolean));
memcpy(data->simulationInfo->relations, optData->re, nRelations*sizeof(modelica_boolean));
memcpy(data->simulationInfo->storedRelations, optData->storeR, nRelations*sizeof(modelica_boolean));

copy_initial_values(optData, data);

for(i = 0, shift = 0; i < nsi-1; ++i){
for(j = 0; j < np; ++j, shift += nv){
Expand Down Expand Up @@ -763,7 +757,6 @@ static inline void updateDOSystem(OptData * optData, DATA * data, threadData_t *
MMC_TRY_INTERNAL(simulationJumpBuffer)
#endif
data->callback->input_function(data, optData->threadData);
/*data->callback->functionDAE(data);*/
updateDiscreteSystem(data, optData->threadData);

if(index){
Expand All @@ -779,7 +772,7 @@ static inline void updateDOSystem(OptData * optData, DATA * data, threadData_t *
* helper optData2ModelData
* author: Vitalij Ruge
**/
static inline void setLocalVars(OptData * optData, DATA * data, const double * const vopt,
void setLocalVars(OptData * optData, DATA * data, const double * const vopt,
const int i, const int j, const int shift){
short l;
int k;
Expand All @@ -790,10 +783,13 @@ static inline void setLocalVars(OptData * optData, DATA * data, const double * c
const int nx = optData->dim.nx;
const int nv = optData->dim.nv;


memcpy(optData->v[i][j], data->localData[l]->realVars, optData->dim.nReal*sizeof(modelica_real));

This comment has been minimized.

Copy link
@mahge

mahge Nov 2, 2021

Contributor

@vruge I think you missed initializing l here. Unfortunately, I do not think you can do the memset() here. Entries of data->localData[:]->realVars are not contiguous :(.

for(l = 0; l < 3; ++l){
data->localData[l]->realVars = optData->v[i][j];
data->localData[l]->timeValue = (modelica_real) optData->time.t[i][j];
}

for(l = 0; l < 2; ++l)
if(optData->s.matrix[l])
data->simulationInfo->analyticJacobians[indexBC[l]].tmpVars = dim->analyticJacobians_tmpVars[l][i][j];
Expand Down Expand Up @@ -952,13 +948,12 @@ static inline void pickUpStates(OptData* optData){
int i, j;
double start_value;
char buffer[200];
const int nReal = optData->data->modelData->nVariablesReal;
rewind(pFile);
for(i =0; i< n; ++i){
fscanf(pFile, "%s", buffer);
if (fscanf(pFile, "%lf", &start_value) <= 0) continue;

for(j = 0, b = 0; j < nReal; ++j){
for(j = 0, b = 0; j < optData->dim.nReal; ++j){
if(!strcmp(optData->data->modelData->realVarsData[j].info.name, buffer)){
optData->data->localData[0]->realVars[j] = start_value;
optData->data->localData[1]->realVars[j] = start_value;
Expand Down
Expand Up @@ -48,6 +48,8 @@ void diffSynColoredOptimizerSystem(OptData *optData, modelica_real **J, const in
void diffSynColoredOptimizerSystemF(OptData *optData, modelica_real **J);
void debugeJac(OptData * optData,Number* vopt);
void debugeSteps(OptData * optData, modelica_real*vopt, modelica_real * lambda);
void copy_initial_values(OptData * optData, DATA* data);
void setLocalVars(OptData * optData, DATA * data, const double * const vopt, const int i, const int j, const int shift);

/*ipopt*/

Expand Down
43 changes: 1 addition & 42 deletions OMCompiler/SimulationRuntime/c/optimization/eval_all/EvalL.c
Expand Up @@ -88,36 +88,6 @@ Bool ipopt_h(int n, double *vopt, Bool new_x, double obj_factor, int m, double *
}
/*******************/
}

#if 0
{
printf("\nnele_hess = %i, %i",nele_hess,k);
FILE *pFile;
char buffer[4096];
pFile = fopen("hesse_struct.m", "wt");
if(pFile == NULL)
printf("\n\nError");
fprintf(pFile, "%s", "clear H\n");
fprintf(pFile, "%s", "%%%%%%%%%%%%%%%%%%%%%%\n");
fprintf(pFile, "%s", "nz = ");
fprintf(pFile, "%i", nele_hess);
fprintf(pFile, "%s", "\nnumberVars = ");
fprintf(pFile, "%i", optData->dim.NV);
fprintf(pFile, "%s", "\nnumberconstraints = ");
fprintf(pFile, "%i", m);
fprintf(pFile, "%s", "\nNumberOfIntervalls = ");
fprintf(pFile, "%i", nsi);
fprintf(pFile, "\nH=sparse(%i,%i);", optData->dim.NV,optData->dim.NV);
fprintf(pFile, "%s", "%%%%%%%%%%%%%%%%%%%%%%\n");
for(i=0; i< k; ++i){
sprintf(buffer, "H(%i,%i) = 1;\n", iRow[i]+1, iCol[i]+1);
fprintf(pFile,"%s", buffer);
}
fprintf(pFile, "%s", "%%%%%%%%%%%%%%%%%%%%%%\n");
fprintf(pFile, "%s", "spy(H)\n");
}
assert(0);
#endif
}else if(keepH){
memcpy(values,optData->oldH,nele_hess*sizeof(double));
}else{
Expand All @@ -127,10 +97,6 @@ Bool ipopt_h(int n, double *vopt, Bool new_x, double obj_factor, int m, double *
const int nJ = optData->dim.nJ;
modelica_boolean upC;
modelica_boolean upC2;
const int nBoolean = optData->data->modelData->nVariablesBoolean;
const int nInteger = optData->data->modelData->nVariablesInteger;
const int nReal = optData->dim.nReal;
const int nRelations = optData->data->modelData->nRelations;
DATA * data = optData->data;

upC = obj_factor != 0;
Expand All @@ -145,14 +111,7 @@ Bool ipopt_h(int n, double *vopt, Bool new_x, double obj_factor, int m, double *
optData->dim.iter_updateHessian = 0;
upC2 = upC && optData->s.mayer;
upC = upC && optData->s.lagrange;

memcpy(data->localData[0]->integerVars, optData->i0, nInteger*sizeof(modelica_integer));
memcpy(data->localData[0]->booleanVars, optData->b0, nBoolean*sizeof(modelica_boolean));
memcpy(data->simulationInfo->integerVarsPre, optData->i0Pre, nInteger*sizeof(modelica_integer));
memcpy(data->simulationInfo->booleanVarsPre, optData->b0Pre, nBoolean*sizeof(modelica_boolean));
memcpy(data->simulationInfo->realVarsPre, optData->v0Pre, nReal*sizeof(modelica_real));
memcpy(data->simulationInfo->relationsPre, optData->rePre, nRelations*sizeof(modelica_boolean));
memcpy(data->simulationInfo->relations, optData->re, nRelations*sizeof(modelica_boolean));
copy_initial_values(optData, data);

for(ii = 0, k = 0, v = vopt, la = lambda; ii + 1 < nsi; ++ii){
for(p = 1; p < np1; ++p, v += nv, la += nJ){
Expand Down
101 changes: 101 additions & 0 deletions testsuite/openmodelica/cruntime/optimization/basic/InputOptIssues.mo
@@ -0,0 +1,101 @@
package InputOptIssues
model Trapezoid "Dynamical Optimization of Ideal Drive"
parameter Real p = 1 "required for optimization";
parameter Real powLim = 9000;
Real power = torque * angSpeed;
//
/*** Optimization requests ***/
input Real torque(min = -90, max = 90);
Real targetPhi(nominal = 1) = -torque2.flange.phi "minimize -pos(tf)" annotation(
isMayer = true);
Real angSpeed(min = 0, max = 0) = p * der(torque2.flange.phi) annotation(
isFinalConstraint = true);
Real pow(min = -powLim, max = powLim) = power annotation(
isConstraint = true);
/*** end of Optimization requests ***/
//
Real constPhi(nominal = 100) = -torque2.flange.phi "minimize -phi(tf)" annotation(
isLagrange = true);
Modelica.Blocks.Sources.RealExpression realexp(y = torque) annotation(
Placement(visible = true, transformation(origin = {0, -34}, extent = {{10, -10}, {-10, 10}}, rotation = 0)));
Modelica.Mechanics.Rotational.Components.Inertia inertia1 annotation(
Placement(visible = true, transformation(origin = {-10, 0}, extent = {{-10, -10}, {10, 10}}, rotation = 0)));
Modelica.Mechanics.Rotational.Sources.Torque torque1 annotation(
Placement(visible = true, transformation(origin = {-48, 0}, extent = {{-10, -10}, {10, 10}}, rotation = 0)));
Modelica.Blocks.Sources.Trapezoid trapezoid1(amplitude = 10, falling = 1, period = 5, rising = 1, startTime = 1, width = 1) annotation(
Placement(visible = true, transformation(origin = {58, 0}, extent = {{10, -10}, {-10, 10}}, rotation = 0)));
Modelica.Mechanics.Rotational.Sources.Torque torque2 annotation(
Placement(visible = true, transformation(origin = {24, 0}, extent = {{10, -10}, {-10, 10}}, rotation = 0)));
equation
connect(trapezoid1.y, torque2.tau) annotation(
Line(points = {{47, 0}, {42, 0}, {42, -2}, {37, -2}}, color = {0, 0, 127}));
connect(torque1.flange, inertia1.flange_a) annotation(
Line(points = {{-38, 0}, {-18, 0}}));
connect(inertia1.flange_b, torque2.flange) annotation(
Line(points = {{2, 0}, {14, 0}}));
connect(realexp.y, torque1.tau) annotation(
Line(points = {{-9, -34}, {-70, -34}, {-70, 0}, {-60, 0}}, color = {0, 0, 127}));
annotation(
Diagram(coordinateSystem(extent = {{-80, -60}, {80, 40}}, preserveAspectRatio = false, initialScale = 0.1, grid = {2, 2})),
Icon(coordinateSystem(extent = {{-80, -60}, {80, 40}}, preserveAspectRatio = false, initialScale = 0.1, grid = {2, 2})),
experiment(StartTime = 0, StopTime = 8, Tolerance = 1e-07, Interval = 0.16),
__OpenModelica_commandLineOptions = "+gDynOpt",
__OpenModelica_simulationFlags(optimizerNP = "1", s = "optimization"));
end Trapezoid;

model TimeTable "Dynamical Optimization of Ideal Drive"
parameter Real p = 1 "required for optimization";
parameter Real powLim = 9000;
Real power = torque * angSpeed;
//
/*** Optimization requests ***/
input Real torque(min = -90, max = 90);
Real targetPhi(nominal = 1) = -torque2.flange.phi "minimize -pos(tf)" annotation(
isMayer = true);
Real angSpeed(min = 0, max = 0) = p * der(torque2.flange.phi) annotation(
isFinalConstraint = true);
Real pow(min = -powLim, max = powLim) = power annotation(
isConstraint = true);
/*** end of Optimization requests ***/
//
Real constPhi(nominal = 100) = -torque2.flange.phi "minimize -phi(tf)" annotation(
isLagrange = true);
Modelica.Blocks.Sources.RealExpression realexp(y = torque) annotation(
Placement(visible = true, transformation(origin = {0, -34}, extent = {{10, -10}, {-10, 10}}, rotation = 0)));
Modelica.Mechanics.Rotational.Components.Inertia inertia1 annotation(
Placement(visible = true, transformation(origin = {-10, 0}, extent = {{-10, -10}, {10, 10}}, rotation = 0)));
Modelica.Mechanics.Rotational.Sources.Torque torque1 annotation(
Placement(visible = true, transformation(origin = {-48, 0}, extent = {{-10, -10}, {10, 10}}, rotation = 0)));
Modelica.Mechanics.Rotational.Sources.Torque torque2 annotation(
Placement(visible = true, transformation(origin = {24, 0}, extent = {{10, -10}, {-10, 10}}, rotation = 0)));
Modelica.Blocks.Sources.TimeTable timeTable1(table = [
0, 0; 1, 0;
2,10; 3, 10;
4, 0; 5, 0;
6,10; 7,10;
8,0; 9,0]
) annotation(
Placement(visible = true, transformation(origin = {64, 0}, extent = {{10, -10}, {-10, 10}}, rotation = 0)));
equation
connect(timeTable1.y, torque2.tau) annotation(
Line(points = {{52, 0}, {38, 0}, {38, 0}, {36, 0}}, color = {0, 0, 127}));
connect(torque1.flange, inertia1.flange_a) annotation(
Line(points = {{-38, 0}, {-18, 0}}));
connect(inertia1.flange_b, torque2.flange) annotation(
Line(points = {{2, 0}, {14, 0}}));
connect(realexp.y, torque1.tau) annotation(
Line(points = {{-9, -34}, {-70, -34}, {-70, 0}, {-60, 0}}, color = {0, 0, 127}));
annotation(
Diagram(coordinateSystem(extent = {{-80, -60}, {80, 40}}, preserveAspectRatio = false, initialScale = 0.1, grid = {2, 2})),
Icon(coordinateSystem(extent = {{-80, -60}, {80, 40}}, preserveAspectRatio = false, initialScale = 0.1, grid = {2, 2})),
experiment(StartTime = 0, StopTime = 10, Tolerance = 1e-07, Interval = 0.2),
__OpenModelica_commandLineOptions = "+gDynOpt",
__OpenModelica_simulationFlags(optimizerNP = "1", s = "optimization"));
end TimeTable;



annotation(
Diagram(coordinateSystem(extent = {{-100, -80}, {100, 80}})),
uses(Modelica(version = "3.2.3")));
end InputOptIssues;

0 comments on commit 998ecd7

Please sign in to comment.