Skip to content

Commit

Permalink
- fixed input update
Browse files Browse the repository at this point in the history 
- improve interpolation for input
- improve initial guess 


git-svn-id: https://openmodelica.org/svn/OpenModelica/trunk@19316 f25d12d1-65f4-0310-ae8a-bbce733d8d8e
  • Loading branch information
Vitalij Ruge committed Feb 26, 2014
1 parent 05b9ab5 commit cb0c363
Show file tree
Hide file tree
Showing 8 changed files with 124 additions and 41 deletions.
68 changes: 41 additions & 27 deletions SimulationRuntime/c/optimization/initialOptimizer/initial_guess.c
View file
Expand Up @@ -141,33 +141,39 @@ static int initial_guess_ipopt_sim(IPOPT_DATA_ *iData,SOLVER_INFO* solverInfo)
data->simulationInfo.inputVars[i] = u0[i];
else
externalInputUpdate(data);

if(iData->preSim){
printf("\npre simulation");
printf("\n========================================================");
printf("\nstart pre simulation");
printf("\n--------------------------------------------------------");
printf("\nfrom %g to %g", iData->t0, iData->startTimeOpt );
pre_ipopt_sim(iData, solverInfo);
printf("\npre simulation done!!!");
printf("\nfinished pre simulation");
printf("\n========================================================\n");
}

printGuess = (short)(ACTIVE_STREAM(LOG_INIT) && !ACTIVE_STREAM(LOG_SOLVER));
if(printGuess){
printf("\n****initial guess****");
printf("\n #####done time[%i] = %f",0,iData->time[0]);
printf("\nInitial Guess");
printf("\n========================================================\n");
printf("\ndone: time[%i] = %g",0,iData->time[0]);
}

for(i=0, k=1, v=iData->v + iData->nv; i<iData->nsi; ++i){
for(jj=0; jj<iData->deg; ++jj, ++k){
a = 1.0;
while(iData->data->localData[0]->timeValue <= iData->time[k]){
solverInfo->currentStepSize = a*(iData->time[k] - iData->time[k-1]);
if(data->simulationInfo.external_input.active)
externalInputUpdate(data);
externalInputUpdate(data);
err = dasrt_step(data, solverInfo);
if(err<0)
a *= 0.5;
if(a < 1e-20)
if((double)a < (double)1e-20)
break;
}

if(printGuess)
printf("\n #####done time[%i] = %f", k, iData->time[k]);
printf("\ndone: time[%i] = %g", k, iData->time[k]);

for(j=0; j< iData->nx; ++j)
v[j] = sData->realVars[j] * iData->scalVar[j];
Expand All @@ -192,8 +198,11 @@ static int initial_guess_ipopt_sim(IPOPT_DATA_ *iData,SOLVER_INFO* solverInfo)
}
}

if(printGuess)
printf("\n*****initial guess done*****");
if(printGuess){
printf("\n--------------------------------------------------------");
printf("\nfinished: Initial Guess");
printf("\n========================================================\n");
}

dasrt_deinitial(solverInfo->solverData);
externalInputFree(data);
Expand All @@ -208,42 +217,47 @@ static int initial_guess_ipopt_sim(IPOPT_DATA_ *iData,SOLVER_INFO* solverInfo)

static int pre_ipopt_sim(IPOPT_DATA_ *iData,SOLVER_INFO* solverInfo)
{
int k = 1,i,j;
int k = 1,i,j,err;
double t;
long double a;
DATA * data = iData->data;

if(iData->time[0] > iData->startTimeOpt)
iData->time[0] = iData->startTimeOpt;

while(iData->data->localData[0]->timeValue < iData->startTimeOpt){
while(iData->data->localData[0]->timeValue <= iData->time[k]){
a = 1.0;
while(iData->data->localData[0]->timeValue < iData->time[k]){
t = iData->time[k];
if(t > iData->startTimeOpt)
t = iData->startTimeOpt;
solverInfo->currentStepSize = t - iData->time[k-1];
solverInfo->currentStepSize = a*(t - iData->time[k-1]);
iData->data->localData[1]->timeValue = iData->time[k-1];
iData->data->localData[0]->timeValue = t;
externalInputUpdate(data);
dasrt_step(data, solverInfo);
err = dasrt_step(data, solverInfo);
if(err<0)
a *= 0.5;
if((double)a < (double)1e-20)
break;
}
data->simulationInfo.terminal = 1;
sim_result.emit(&sim_result,data);
data->simulationInfo.terminal = 0;
rotateRingBuffer(iData->data->simulationData, 1, (void**) iData->data->localData);
++k;
}

if(iData->data->localData[0]->timeValue > iData->startTimeOpt){
solverInfo->currentStepSize = iData->startTimeOpt - iData->time[k-1];
iData->data->localData[1]->timeValue = iData->time[k-1];
externalInputUpdate(data);
dasrt_step(data, solverInfo);
data->simulationInfo.terminal = 1;
sim_result.emit(&sim_result,data);
data->simulationInfo.terminal = 0;
rotateRingBuffer(iData->data->simulationData, 1, (void**) iData->data->localData);
}

iData->t0 = iData->data->localData[0]->timeValue;
printf("time = %g | %g", iData->t0,iData->startTimeOpt );
/*ToDo*/
for(i=0; i< iData->nx; ++i)
{
Expand Down Expand Up @@ -278,18 +292,18 @@ static int optimizer_time_setings_update(IPOPT_DATA_ *iData)
for(i = 0,k=0,id=0; i<iData->nsi; ++i,id += iData->deg){
if(i){
if(iData->deg == 3){
iData->time[++k] = iData->time[id] + iData->c1*iData->dt[i];
iData->time[++k] = iData->time[id] + iData->c2*iData->dt[i];
}
iData->time[++k] = iData->time[id] + iData->c1*iData->dt[i];
iData->time[++k] = iData->time[id] + iData->c2*iData->dt[i];
}
iData->time[++k] = iData->time[0]+ (i+1)*iData->dt[i];
}else{
if(iData->deg == 3){
iData->time[++k] = iData->time[id] + iData->e1*iData->dt[i];
iData->time[++k] = iData->time[id] + iData->e2*iData->dt[i];
}
}else{
if(iData->deg == 3){
iData->time[++k] = iData->time[id] + iData->e1*iData->dt[i];
iData->time[++k] = iData->time[id] + iData->e2*iData->dt[i];
}
iData->time[++k] = iData->time[0]+ (i+1)*iData->dt[i];
}
}
}
}
}
iData->time[k] = iData->tf;
return 0;
Expand Down
5 changes: 3 additions & 2 deletions SimulationRuntime/c/optimization/mainOptimizer/ipoptODE.c
View file
Expand Up @@ -209,8 +209,8 @@ static int res2file(IPOPT_DATA_ *iData,SOLVER_INFO* solverInfo)

fprintf(pFile, "%lf ",iData->time[iData->current_time]);
for(i=0,j=iData->nx; i< iData->nu; ++i,++j){
data->simulationInfo.inputVars[i] = iData->v[k++]*iData->vnom[j];
fprintf(pFile, "%lf ", data->simulationInfo.inputVars[i]);
data->simulationInfo.inputVars[i] = iData->v[k]*iData->vnom[j];
fprintf(pFile, "%lf ", iData->v[k++]*iData->vnom[j]);
}
fprintf(pFile, "%s", "\n");

Expand Down Expand Up @@ -298,6 +298,7 @@ static int set_optimizer_flags(IPOPT_DATA_ *iData, IpoptProblem *nlp)
/*
* AddIpoptStrOption(nlp, "derivative_test_print_all", "yes");
* AddIpoptNumOption(nlp,"derivative_test_perturbation",1e-6);
*
*/

AddIpoptIntOption(*nlp, "max_iter", 5000);
Expand Down
28 changes: 25 additions & 3 deletions SimulationRuntime/c/simulation/solver/dassl.c
View file
Expand Up @@ -308,8 +308,7 @@ int dasrt_step(DATA* simData, SOLVER_INFO* solverInfo)

/* read input vars */
if(solverInfo->solverMethod != S_OPTIMIZATION) {
if(simData->simulationInfo.external_input.active)
externalInputUpdate(simData);
externalInputUpdate(simData);
simData->callback->input_function(simData);
}

Expand Down Expand Up @@ -374,6 +373,11 @@ int dasrt_step(DATA* simData, SOLVER_INFO* solverInfo)
fflush(stderr);
fflush(stdout);
retVal = continue_DASRT(&dasslData->idid, &simData->simulationInfo.tolerance);
/* read input vars */
if(solverInfo->solverMethod != S_OPTIMIZATION) {
externalInputUpdate(simData);
simData->callback->input_function(simData);
}
simData->callback->functionODE(simData);
warningStreamPrint(LOG_STDOUT, 0, "can't continue. time = %f", sData->timeValue);
return retVal;
Expand Down Expand Up @@ -495,7 +499,13 @@ int functionODE_residual(double *t, double *y, double *yd, double *delta,
#if !defined(OMC_EMCC)
MMC_TRY_INTERNAL(simulationJumpBuffer)
#endif
data->callback->functionODE(data);

/* read input vars */
externalInputUpdate(data);
data->callback->input_function(data);

/* eval input vars */
data->callback->functionODE(data);

/* get the difference between the temp_xd(=localData->statesDerivatives)
and xd(=statesDerivativesBackup) */
Expand Down Expand Up @@ -532,6 +542,10 @@ int function_ZeroCrossingsDASSL(fortran_integer *neqm, double *t, double *y,
timeBackup = data->localData[0]->timeValue;

data->localData[0]->timeValue = *t;
/* read input vars */
externalInputUpdate(data);
data->callback->input_function(data);
/* eval ode*/
data->callback->functionODE(data);
data->callback->functionAlgebraics(data);

Expand Down Expand Up @@ -665,6 +679,10 @@ static int JacobianSymbolicColored(double *t, double *y, double *yprime, double

data->localData[0]->timeValue = *t;
data->localData[0]->realVars = y;
/* read input vars */
externalInputUpdate(data);
data->callback->input_function(data);
/* eval ode*/
data->callback->functionODE(data);
functionJacAColored(data, pd);

Expand Down Expand Up @@ -698,6 +716,10 @@ static int JacobianSymbolic(double *t, double *y, double *yprime, double *deltaD

data->localData[0]->timeValue = *t;
data->localData[0]->realVars = y;
/* read input vars */
externalInputUpdate(data);
data->callback->input_function(data);
/* eval ode*/
data->callback->functionODE(data);
functionJacASym(data, pd);

Expand Down
4 changes: 4 additions & 0 deletions SimulationRuntime/c/simulation/solver/events.c
View file
Expand Up @@ -428,6 +428,10 @@ double bisection(DATA* data, double* a, double* b, double* states_a, double* sta
}

/*calculates Values dependents on new states*/
/* read input vars */
externalInputUpdate(data);
data->callback->input_function(data);
/* eval ODE*/
data->callback->functionODE(data);
data->callback->functionAlgebraics(data);

Expand Down
45 changes: 42 additions & 3 deletions SimulationRuntime/c/simulation/solver/external_input.c
View file
Expand Up @@ -117,24 +117,63 @@ int externalInputFree(DATA* data)
int externalInputUpdate(DATA* data)
{
double u1, u2;
double t, t1, dt;
double t, t1, t2, dt;
int i;

if(!data->simulationInfo.external_input.active)
return -1;

t = data->localData[0]->timeValue;
t1 = data->simulationInfo.external_input.t[data->simulationInfo.external_input.i];
t2 = data->simulationInfo.external_input.t[data->simulationInfo.external_input.i+1];

if(t == t1){
for(i = 0; i < data->modelData.nInputVars; ++i){
data->simulationInfo.inputVars[i] = data->simulationInfo.external_input.u[data->simulationInfo.external_input.i][i];
}
return 1;
}else if(t == t2){
for(i = 0; i < data->modelData.nInputVars; ++i){
data->simulationInfo.inputVars[i] = data->simulationInfo.external_input.u[data->simulationInfo.external_input.i][i+1];
}
return 1;
}

while(i > 0 && t < t1){
--data->simulationInfo.external_input.i;
t1 = data->simulationInfo.external_input.t[data->simulationInfo.external_input.i];
}

while(t > data->simulationInfo.external_input.t[data->simulationInfo.external_input.i+1]
&& data->simulationInfo.external_input.i+1 < (data->simulationInfo.external_input.n-1)){
++data->simulationInfo.external_input.i;
}

if(t == t1){
for(i = 0; i < data->modelData.nInputVars; ++i){
data->simulationInfo.inputVars[i] = data->simulationInfo.external_input.u[data->simulationInfo.external_input.i][i];
}
return 1;
}else if(t == t2){
for(i = 0; i < data->modelData.nInputVars; ++i){
data->simulationInfo.inputVars[i] = data->simulationInfo.external_input.u[data->simulationInfo.external_input.i][i+1];
}
return 1;
}

dt = (data->simulationInfo.external_input.t[data->simulationInfo.external_input.i+1] - data->simulationInfo.external_input.t[data->simulationInfo.external_input.i]);
t1 = data->simulationInfo.external_input.t[data->simulationInfo.external_input.i];
for(i = 0; i < data->modelData.nInputVars; ++i){
u1 = data->simulationInfo.external_input.u[data->simulationInfo.external_input.i][i];
u2 = data->simulationInfo.external_input.u[data->simulationInfo.external_input.i+1][i];
data->simulationInfo.inputVars[i] = (u1*(dt+t1-t)+(t-t1)*u2)/dt;

if(u1 != u2){
if(sign(u1) == sign(u2))
data->simulationInfo.inputVars[i] = ((u1*(dt+t1) + t*u2) - (t*u1+t1*u2))/dt;
else
data->simulationInfo.inputVars[i] = (u1*(dt+t1-t)+(t-t1)*u2)/dt;
}else{
data->simulationInfo.inputVars[i] = u1;
}
}
return 0;
}
1 change: 1 addition & 0 deletions SimulationRuntime/c/simulation/solver/perform_simulation.c
View file
Expand Up @@ -58,6 +58,7 @@
*/
void updateContinuousSystem(DATA *data)
{
externalInputUpdate(data);
data->callback->input_function(data);
data->callback->functionODE(data);
data->callback->functionAlgebraics(data);
Expand Down
3 changes: 3 additions & 0 deletions SimulationRuntime/c/simulation/solver/radau.c
View file
Expand Up @@ -443,6 +443,9 @@ static int refreshModell(DATA* data, double* x, double time)

memcpy(sData->realVars, x, sizeof(double)*data->modelData.nStates);
sData->timeValue = time;
/* read input vars */
externalInputUpdate(data);
data->callback->input_function(data);
data->callback->functionODE(data);

return 0;
Expand Down
11 changes: 5 additions & 6 deletions SimulationRuntime/c/simulation/solver/solver_main.c
View file
Expand Up @@ -338,7 +338,12 @@ int initializeModel(DATA* data, const char* init_initMethod,
copyStartValuestoInitValues(data);

/* read input vars */
externalInputUpdate(data);
data->callback->input_function(data);
/* update start values for inputs if input is set */
if(data->simulationInfo.external_input.active){
data->callback->input_function_init(data);
}

data->localData[0]->timeValue = simInfo->startTime;

Expand Down Expand Up @@ -555,12 +560,6 @@ int solver_main(DATA* data, const char* init_initMethod,
/* allocate SolverInfo memory */
retVal = initializeSolverData(data, &solverInfo);
omc_alloc_interface.collect_a_little();

/* set inputs from file */
if(data->simulationInfo.external_input.active){
externalInputUpdate(data);
data->callback->input_function_init(data);
}

/* initialize all parts of the model */
if(0 == retVal) {
Expand Down

0 comments on commit cb0c363

Please sign in to comment.