Skip to content

Commit

Permalink
-imporved error msg for dyn. optimization
Browse files Browse the repository at this point in the history 
-some fixes from static analysis 


git-svn-id: https://openmodelica.org/svn/OpenModelica/trunk@23004 f25d12d1-65f4-0310-ae8a-bbce733d8d8e
  • Loading branch information
Vitalij Ruge committed Oct 28, 2014
1 parent 090e12d commit f2f0c0c
Show file tree
Hide file tree
Showing 5 changed files with 90 additions and 41 deletions.
1 change: 1 addition & 0 deletions SimulationRuntime/c/optimization/OptimizerData.h
View file
Expand Up @@ -183,6 +183,7 @@ typedef struct OptData{
FILE * pFile;

double *oldH;
int iter_;

}OptData;

Expand Down
5 changes: 3 additions & 2 deletions SimulationRuntime/c/optimization/eval_all/EvalF.c
View file
Expand Up @@ -80,7 +80,7 @@ Bool evalfF(Index n, Number * vopt, Bool new_x, Number *objValue, void * useData
}

i = nsi - 1;
for(j = 0; j< np; ++j)
for(j = 0; j < np; ++j)
erg1 += b[j]*v[i][j][il];

lagrange = (erg*dt[0] + erg1*dt[1]);
Expand Down Expand Up @@ -144,8 +144,9 @@ Bool evalfDiffF(Index n, double * vopt, Bool new_x, Number *gradF, void * useDat
modelica_real * gradM = optData->J[nsi - 1][np -1][nJ1];
if(la){
int i;
const int nnv = n - nv;
for(i = 0; i < nv; ++i)
gradF[n-nv + i] += gradM[i];
gradF[nnv + i] += gradM[i];
}else{
memcpy(gradF + n - nv, gradM, nv*sizeof(modelica_real));
}
Expand Down
111 changes: 77 additions & 34 deletions SimulationRuntime/c/optimization/eval_all/EvalG.c
View file
Expand Up @@ -65,6 +65,8 @@ static inline void structJacC(const modelica_real *const J, double *values,
static inline void printMaxError(Number *g, const int m, const int nx, const int nJ, long double **t,
const int np, const int nsi, DATA * data, OptData * optData);

static inline void debugeJac(OptData * optData);


/* eval constraints
* author: Vitalij Ruge
Expand All @@ -87,8 +89,10 @@ Bool evalfG(Index n, double * vopt, Bool new_x, int m, Number *g, void * useData
double * vv[np+1];
int i, j, k, shift;

if(new_x)

if(new_x){
optData2ModelData(optData, vopt, 1);
}

v = optData->v;
memcpy(a ,optData->rk.a, sizeof(optData->rk.a));
Expand Down Expand Up @@ -151,7 +155,6 @@ Bool evalfG(Index n, double * vopt, Bool new_x, int m, Number *g, void * useData
const int nJ = optData->dim.nJ;
printMaxError(g, m, nx, nJ, optData->time.t, np ,nsi ,optData->data, optData);
}

return TRUE;
}

Expand Down Expand Up @@ -204,9 +207,10 @@ Bool evalfDiffG(Index n, double * vopt, Bool new_x, Index m, Index njac, Index *
modelica_boolean ** J = optData->s.JderCon;
modelica_boolean ** Jf = optData->s.J[2];
int i, j, k, l, ii, cindex;

if(new_x)
++optData->iter_;
if(new_x){
optData2ModelData(optData, vopt, 1);
}
if(np == 3){
/*****************************/
for(j = 0, k = 0; j < np; ++j){
Expand Down Expand Up @@ -315,6 +319,9 @@ Bool evalfDiffG(Index n, double * vopt, Bool new_x, Index m, Index njac, Index *
printf("\nvalues[%i] = %g",i,values[i]);
assert(0);
*/
#if 0
debugeJac(optData);
#endif
}

return TRUE;
Expand Down Expand Up @@ -660,10 +667,10 @@ static inline void printMaxError(Number *g, const int m, const int nx, const int
int index_x = 0;
double gmax = -1;
int i, j, k, l;
int ii, jj, kk = -1;
int ii=-1, jj=-1, kk = -1;
double tmp, tmp1;

for(i = 0, l = 0; i < nsi-1; ++i){
for(i = 0, l = 0; i < nsi; ++i){
for(j = 0; j < np; ++j){
for(k=0; k< nx; ++k){
tmp = fabs(g[l++]);
Expand All @@ -675,9 +682,9 @@ static inline void printMaxError(Number *g, const int m, const int nx, const int
}
}
for(; k< nJ; ++k, ++l){
tmp1 = g[l] > optData->ipop.gmax[l] ? fabs(g[l] - optData->ipop.gmax[l]) : 0.0;
tmp = g[l] < optData->ipop.gmin[l] ? fabs(g[l] - optData->ipop.gmin[l]) : 0.0;
tmp = fmax(tmp, tmp1);
tmp1 = g[l] - optData->ipop.gmax[l]; // > 0
tmp = optData->ipop.gmin[l] - g[l]; // >0
tmp = fmaxl(fmaxl(tmp, tmp1),0.0);
if(tmp > gmax){
ii = i;
jj = j;
Expand All @@ -688,37 +695,73 @@ static inline void printMaxError(Number *g, const int m, const int nx, const int
}
}

for(j = 0; j < np; ++j){
for(k=0; k< nx; ++k){
tmp = fabs(g[l++]);
if(tmp > gmax){
ii = nsi - 1;
jj = j;
kk = k;
gmax = tmp;
}
}
/*
for(; k< nJ+optData->dim.ncf; ++k, ++l){
tmp = fmax(fabs(g[l] - optData->ipop.gmax[l]), fabs(g[l] - optData->ipop.gmin[l]));
if(tmp > gmax){
ii = nsi- 1;
jj = j;
kk = k;
gmax = tmp;
}
/*final constraints*/
for(k=nJ; k< nJ+optData->dim.ncf; ++k, ++l){
tmp1 = g[l] - optData->ipop.gmax[l]; // > 0
tmp = optData->ipop.gmin[l] - g[l]; // >0
tmp = fmaxl(fmaxl(tmp, tmp1),0.0);
if(tmp > gmax){
ii = nsi- 1;
jj = np-1;
kk = k;
gmax = tmp;
}
*/
}
if(k>0){

if(kk>-1){
if(kk < nx){
printf("\nmax error for |%s(%g) - collocation_poly| = %g\n",
data->modelData.realVarsData[kk].info.name, (double)t[ii][jj], gmax);
printf("\nmax error is %g for the approximation of the state %s(time = %g)\n",
gmax, data->modelData.realVarsData[kk].info.name, (double)t[ii][jj]);
}else if(kk < nJ){
printf("\nmax error for |cosntrain[%i](%g)| = %g\n", kk - nx, (double)t[ii][jj], gmax);
const int ll = kk - nx + optData->dim.index_con;
printf("\nmax violation is %g for the constraint %s(time = %g)\n",
gmax, data->modelData.realVarsData[ll].info.name, (double)t[ii][jj]);
}else{
printf("\nmax error for |final_cosntrain[%i](%g)| = %g\n", kk - nJ, (double)t[ii][jj], gmax);
const int ll = kk - nx + optData->dim.index_con;
printf("\nmax violation is %g for the final constraint %s(time = %g)\n", gmax, data->modelData.realVarsData[ll].info.name, (double)t[ii][jj]);
}
}
}

/*!
* generated csv and python script for jacobian
* author: Vitalij Ruge
**/
static inline void debugeJac(OptData * optData){
int i,j,k, jj;
const int nv = optData->dim.nv;
const int nx = optData->dim.nx;
const int nu = optData->dim.nu;
const int nsi = optData->dim.nsi;
const int nJ = optData->dim.nJ;
const int np = optData->dim.np;
const int npv = np*nv;
double **J;

FILE *pFile;
char buffer[4096];

sprintf(buffer, "jac_ana_step_%i.csv", optData->iter_);
pFile = fopen(buffer, "wt");

fprintf(pFile,"name;time;");
for(j = 0; j < nx; ++j)
fprintf(pFile,"%s;",optData->data->modelData.realVarsData[j].info.name);
for(j = 0; j < nu; ++j)
fprintf(pFile, "%s;", optData->dim.inputName[j]);
fprintf(pFile,"\n");

for(i=0;i < nsi; ++i){
for(j = 0; j < nJ; ++j){
J = optData->J[i][j];
for(k = 0; k < nx; ++k){
fprintf(pFile,"%s;%g;",optData->data->modelData.realVarsData[k].info.name,(double)optData->time.t[i][j]);
for(jj = 0; jj < nv; ++jj)
fprintf(pFile,"%g;", J[k][jj]);
fprintf(pFile,"\n");
}
}
}

fclose(pFile);
}
2 changes: 1 addition & 1 deletion SimulationRuntime/c/optimization/eval_all/EvalL.c
View file
Expand Up @@ -226,7 +226,7 @@ static inline void num_hessian0(double * v, const double * const lambda,
v_save = (long double) v[ii];
h = (long double)DF_STEP(v_save);
v[ii] += h;
if( v[ii] > vmax[ii]){
if( v[ii] >= vmax[ii]){
h *= -1.0;
v[ii] = v_save + h;
}
Expand Down
12 changes: 8 additions & 4 deletions SimulationRuntime/c/optimization/optimizer_main.c
View file
Expand Up @@ -100,7 +100,8 @@ static inline void optimizationWithIpopt(OptData*optData){

/*tol */
AddIpoptNumOption(nlp, "tol", optData->data->simulationInfo.tolerance);

AddIpoptStrOption(nlp, "evaluate_orig_obj_at_resto_trial", "yes");

/* print level */
if(ACTIVE_STREAM(LOG_IPOPT_FULL)){
AddIpoptIntOption(nlp, "print_level", 7);
Expand Down Expand Up @@ -143,7 +144,8 @@ static inline void optimizationWithIpopt(OptData*optData){
cflags = (char*)omc_flagValue[FLAG_LS_IPOPT];
if(cflags)
AddIpoptStrOption(nlp, "linear_solver", cflags);

AddIpoptNumOption(nlp,"mumps_pivtolmax",1e-5);


/* max iter */
cflags = (char*)omc_flagValue[FLAG_IPOPT_MAX_ITER];
Expand Down Expand Up @@ -207,9 +209,11 @@ static inline void optimizationWithIpopt(OptData*optData){
/********************************************************************/


if(max_iter >=0)
if(max_iter >=0){
optData->iter_ = 0.0;
res = IpoptSolve(nlp, vopt, NULL, &obj, mult_g, mult_x_L, mult_x_U, (void*)optData);
if(res < 0 && !ACTIVE_STREAM(LOG_IPOPT))
}
if(res != 0 && !ACTIVE_STREAM(LOG_IPOPT))
warningStreamPrint(LOG_STDOUT, 0, "No optimal solution found!\nUse -lv=LOG_IPOPT for more information.");
FreeIpoptProblem(nlp);
}
Expand Down

0 comments on commit f2f0c0c

Please sign in to comment.