Skip to content

Commit

Permalink
- independent memory for jacobian in optimization
Browse files Browse the repository at this point in the history
git-svn-id: https://openmodelica.org/svn/OpenModelica/trunk@19894 f25d12d1-65f4-0310-ae8a-bbce733d8d8e
  • Loading branch information
Vitalij Ruge committed Apr 1, 2014
1 parent 1b0cfe6 commit 5b56d62
Show file tree
Hide file tree
Showing 4 changed files with 42 additions and 27 deletions.
22 changes: 11 additions & 11 deletions SimulationRuntime/c/optimization/constraints/evalfDiffG.c
Expand Up @@ -119,50 +119,50 @@ Bool evalfDiffG(Index n, double * x, Bool new_x, Index m, Index njac, Index *iRo
tmp[1] = dtime->dt[0]*mbase->d[1][4];
tmp[2] = dtime->dt[0]*mbase->d[2][4];

diff_functionODE(x, dtime->t0 , iData, df->J0);
diff_functionODE(x, dtime->t0 , iData, df->J[0]);

for(i = 0, id = dim->nv, k = 0; i<1; ++i){
tmp_index = i*dim->deg;
for(l=0; l<dim->deg; ++l, id += dim->nv){
diff_functionODE(x+id , dtime->time[tmp_index + l] , iData, iData->df.J);
diff_functionODE(x+id , dtime->time[tmp_index + l] , iData, iData->df.J[i+1]);
for(j=0; j<dim->nx; ++j){
switch(l){
case 0:
lobattoJac1(mbase->d[l], df->J[j], df->J0[j], dtime->dt[i], values, dim->nv, &k, j, tmp[l], iData);
lobattoJac1(mbase->d[l], df->J[i+1][j], df->J[0][j], dtime->dt[i], values, dim->nv, &k, j, tmp[l], iData);
break;
case 1:
lobattoJac2(mbase->d[l], df->J[j], df->J0[j], dtime->dt[i], values, dim->nv, &k, j, tmp[l], iData);
lobattoJac2(mbase->d[l], df->J[i+1][j], df->J[0][j], dtime->dt[i], values, dim->nv, &k, j, tmp[l], iData);
break;
case 2:
lobattoJac3(mbase->d[l], df->J[j], df->J0[j], dtime->dt[i], values, dim->nv, &k, j, tmp[l], iData);
lobattoJac3(mbase->d[l], df->J[i+1][j], df->J[0][j], dtime->dt[i], values, dim->nv, &k, j, tmp[l], iData);
break;
}
}
for(;j<nng; ++j){
conJac(df->J[j], values, dim->nv, &k, j, iData);
conJac(df->J[i][j], values, dim->nv, &k, j, iData);
}
}
}

for(; i<dim->nsi; ++i){
tmp_index = i*iData->dim.deg;
for(l=0; l<dim->deg; ++l, id += dim->nv){
diff_functionODE(x+id, dtime->time[tmp_index + l], iData, df->J);
diff_functionODE(x+id, dtime->time[tmp_index + l], iData, df->J[i+1]);
for(j=0; j<dim->nx; ++j){
switch(l){
case 0:
radauJac1(mbase->a[l], df->J[j], dtime->dt[i], values, dim->nv, &k, j, iData);
radauJac1(mbase->a[l], df->J[i+1][j], dtime->dt[i], values, dim->nv, &k, j, iData);
break;
case 1:
radauJac2(mbase->a[l], df->J[j], dtime->dt[i], values, dim->nv, &k, j, iData);
radauJac2(mbase->a[l], df->J[i+1][j], dtime->dt[i], values, dim->nv, &k, j, iData);
break;
case 2:
radauJac3(mbase->a[l], df->J[j], dtime->dt[i], values, dim->nv, &k, j, iData);
radauJac3(mbase->a[l], df->J[i+1][j], dtime->dt[i], values, dim->nv, &k, j, iData);
break;
}
}
for(;j<nng; ++j)
conJac(df->J[j], values, dim->nv, &k, j, iData);
conJac(df->J[i][j], values, dim->nv, &k, j, iData);
}
}
/*assert(k == njac);*/
Expand Down
Expand Up @@ -101,20 +101,26 @@ int allocateIpoptData(IPOPT_DATA_ *iData)
iData->dtime.time = (long double*)malloc((dim->deg*dim->nsi +1) *sizeof(long double));
iData->start_u = (double*)malloc(dim->nv*sizeof(double));

df->J = (long double**) malloc(dim->nJ * sizeof(long double*));
for(i = 0; i < dim->nJ; i++)
df->J[i] = (long double*) calloc(dim->nv, sizeof(long double));
df->J = (long double***) malloc((dim->nsi*dim->deg + 1) * sizeof(long double**));
for(j = 0; j < (dim->nsi*dim->deg + 1); ++j){
df->J[j] = (long double**) malloc(dim->nJ * sizeof(long double*));
for(i = 0; i < dim->nJ; i++)
df->J[j][i] = (long double*) calloc(dim->nv, sizeof(long double));
}

df->Jh = (long double***) malloc(2 * sizeof(long double**));
for(j = 0; j < 2; ++j){
df->Jh[j] = (long double**) malloc(dim->nJ * sizeof(long double*));
for(i = 0; i < dim->nJ; i++)
df->Jh[j][i] = (long double*) calloc(dim->nv, sizeof(long double));
}

for(i = 0; i< 4 ; ++i)
df->gradF[i] = (long double*) calloc(dim->nv, sizeof(long double));

iData->sv = (double*)malloc(dim->nv*sizeof(double));
iData->sh = (double*)malloc(dim->nJ*sizeof(double));

df->J0 = (long double**) malloc(dim->nJ * sizeof(long double*));
for(i = 0; i < dim->nJ; i++)
df->J0[i] = (long double*) calloc(dim->nv, sizeof(long double));

iData->scaling.scaldt = (long double**) malloc(dim->nx * sizeof(long double*));
for(i = 0; i < dim->nx; i++)
iData->scaling.scaldt[i] = (long double*) calloc(dim->nsi, sizeof(long double));
Expand Down Expand Up @@ -178,11 +184,21 @@ static int freeIpoptData(IPOPT_DATA_ *iData)
OPTIMIZER_BOUNDS *bounds = &iData->bounds;
OPTIMIZER_DF *df = &iData->df;
for(i = 0; i < dim->nJ; i++){
free(df->J[i]);
free(df->J0[i]);
free(sopt->knowedJ[i]);
}
free(df->J0);

for(j = 0; j < (dim->nsi*dim->deg + 1); ++j){
for(i = 0; i < dim->nJ; i++)
free(df->J[j][i]);
free(df->J[j]);
}

for(j = 0; j < 2; ++j){
for(i = 0; i < dim->nJ; i++)
free(df->Jh[j][i]);
free(df->Jh[j]);
}

free(df->J);
free(sopt->knowedJ);

Expand Down
5 changes: 2 additions & 3 deletions SimulationRuntime/c/optimization/ipoptODEstruct.h
Expand Up @@ -167,9 +167,8 @@ typedef struct OPTIMIZER_SCALING{

typedef struct OPTIMIZER_DF{

long double **J0;
long double **J;

long double ***J;
long double ***Jh;
long double ** gradFomc;

long double * gradF[4];
Expand Down
6 changes: 3 additions & 3 deletions SimulationRuntime/c/optimization/lagrangeFun/ipopt_hessian.c
Expand Up @@ -251,7 +251,7 @@ static int num_hessian(double *v, double t, IPOPT_DATA_ *iData, double *lambda,
OPTIMIZER_DIM_VARS *dim = &iData->dim;
int nJ = (t>(double)iData->dtime.t0) ? dim->nx + dim->nc : dim->nx;

diff_functionODE(v, t , iData, iData->df.J0);
diff_functionODE(v, t , iData, iData->df.Jh[0]);
upCost = (lagrange_yes || mayer_yes) && (obj_factor!=0);

if(upCost)
Expand All @@ -261,7 +261,7 @@ static int num_hessian(double *v, double t, IPOPT_DATA_ *iData, double *lambda,
v_save = (long double)v[i];
h = (long double)DF_STEP(v_save, iData->scaling.vnom[i]);
v[i] += h;
diff_functionODE(v, t , iData, iData->df.J);
diff_functionODE(v, t , iData, iData->df.Jh[1]);

if(upCost)
updateCost(v,t,iData,lagrange_yes,mayer_yes, iData->df.gradF[0], iData->df.gradF[1]);
Expand All @@ -272,7 +272,7 @@ static int num_hessian(double *v, double t, IPOPT_DATA_ *iData, double *lambda,
if(iData->sopt.Hg[i][j]){
for(l = 0; l< nJ; ++l){
if(iData->sopt.knowedJ[l][j] + iData->sopt.knowedJ[l][i] >= 2 && lambda[l] != 0.0)
iData->df.H[l][i][j] = (long double)(iData->df.J[l][j] - iData->df.J0[l][j])*lambda[l]/h;
iData->df.H[l][i][j] = (long double)(iData->df.Jh[1][l][j] - iData->df.Jh[0][l][j])*lambda[l]/h;
else
iData->df.H[l][i][j] = (long double) 0.0;
iData->df.H[l][j][i] = iData->df.H[l][i][j];
Expand Down

0 comments on commit 5b56d62

Please sign in to comment.