Skip to content

Commit

Permalink
- added dense numerical jac
Browse files Browse the repository at this point in the history 
git-svn-id: https://openmodelica.org/svn/OpenModelica/trunk@19089 f25d12d1-65f4-0310-ae8a-bbce733d8d8e
  • Loading branch information
Vitalij Ruge committed Feb 13, 2014
1 parent b6f1f20 commit 05c9498
Show file tree
Hide file tree
Showing 4 changed files with 79 additions and 34 deletions.
27 changes: 13 additions & 14 deletions SimulationRuntime/c/optimization/constraints/evalfG.c
View file
Expand Up @@ -138,28 +138,27 @@ int functionODE_(double * x, double *u, double t, double * dotx, IPOPT_DATA_ *iD
**/
int diff_functionODE(double* v, double t, IPOPT_DATA_ *iData, double **J)
{
double h;
int i, j, k;
double vsave;
double tmp;
int i, j;
double *x, *u;
long double rcal;
int nJ = iData->nx + iData->nc;
int nJ = (int)iData->nJ;
x = v;
u = v + iData->nx;

if(iData->useNumJac){
num_diff_symColoredODE(v,t,iData,J);
for(i = 0;i<iData->nv;++i)

num_diff_symColoredODE(v,t,iData,J);
for(i = 0;i<iData->nv;++i)
for(j = 0; j <iData->nx; ++j)
iData->numJ[j][i] *= iData->scalf[j];

}else{
refreshSimData(x,u,t,iData);
refreshSimData(x,u,t,iData);
diff_symColoredODE(v,t,iData,J);
for(i = 0;i<iData->nv;++i){
for(j = 0; j <iData->nx; ++j)
J[j][i] *= iData->scalf[j]*iData->vnom[i];
for(; j <nJ; ++j)
J[j][i] *= iData->vnom[i];
for(j = 0; j <iData->nx; ++j)
J[j][i] *= iData->scalf[j]*iData->vnom[i];
for(; j <nJ; ++j)
J[j][i] *= iData->vnom[i];
}
}

Expand Down Expand Up @@ -202,7 +201,7 @@ int num_diff_symColoredODE(double *v, double t, IPOPT_DATA_ *iData, double **J)
memcpy(iData->vsave, v, sizeof(double)*nx);

for(ii = 0; ii<nx; ++ii){
iData->eps[ii] = DF_STEP(v[ii], iData->vnom[ii]);;
iData->eps[ii] = DF_STEP(v[ii], iData->vnom[ii]);
}

for(i = 1; i < data->simulationInfo.analyticJacobians[index].sparsePattern.maxColors + 1; ++i){
Expand Down
73 changes: 65 additions & 8 deletions SimulationRuntime/c/optimization/initialOptimizer/allocate_ipopt.c
View file
Expand Up @@ -38,13 +38,15 @@

#include "../ipoptODEstruct.h"
#include "simulation_data.h"
#include "../../simulation/options.h"

#ifdef WITH_IPOPT

#include "../localFunction.h"

static int local_jac_struct(IPOPT_DATA_ *iData, int *nng);
static int local_jac_struct_print(IPOPT_DATA_ *iData);
static int local_jac_struct_dense(IPOPT_DATA_ *iData, int * nng);
static int check_nominal(IPOPT_DATA_ *iData, double min, double max, double nominal, short set, int i, double x0);
static int optimizer_coeff_setings(IPOPT_DATA_ *iData);
static int optimizer_bounds_setings(DATA *data, IPOPT_DATA_ *iData);
Expand Down Expand Up @@ -317,7 +319,31 @@ int loadDAEmodel(DATA *data, IPOPT_DATA_ *iData)
optimizer_coeff_setings(iData);

/***********************/
local_jac_struct(iData, &id);
char *cflags;

iData->useNumJac = 0;
cflags = (char*)omc_flagValue[FLAG_IPOPT_JAC];
if(cflags){
if(!strcmp(cflags,"NUM"))
iData->useNumJac = 1;
else if(!strcmp(cflags,"SYM") || !strcmp(cflags,"sym"))
iData->useNumJac = 0;
else if(!strcmp(cflags,"NUMDENSE") || !strcmp(cflags,"NUMDENSE"))
iData->useNumJac = 2;
else
warningStreamPrint(LOG_STDOUT, 1, "not support ipopt_hesse=%s",cflags);
}

printf("\nnumJ = %i",iData->useNumJac);
if(iData->useNumJac == 2){
local_jac_struct_dense(iData, &id);
printf("\nid = %i", id);
} else {
local_jac_struct(iData, &id);
printf("\nsid = %i", id);
}


iData->njac = iData->deg*(iData->nlocalJac-iData->nx+iData->nsi*iData->nlocalJac+iData->deg*iData->nsi*iData->nx)-iData->deg*id;
iData->nhess = 0.5*iData->nv*(iData->nv + 1)*(1+iData->deg*iData->nsi);

Expand Down Expand Up @@ -362,6 +388,37 @@ int move_grid(IPOPT_DATA_ *iData)
return 0;
}

/*
* function calculates a jacobian matrix struct
* author: vitalij
*/
static int local_jac_struct_dense(IPOPT_DATA_ *iData, int * nng)
{
int **J;
short **Hg;
int i,j,ii;
int nJ = (int)(iData->nJ);

J = iData->knowedJ;
Hg = iData->Hg;
//iData->useNumJac == 2

for(i = 0; i < iData->nJ; ++i)
for(j = 0; j < iData->nv; ++j){
J[i][j] = 1.0;
++iData->nlocalJac;
}


for(i = 0; i <iData->nv; ++i)
for(j = 0; j < iData->nv; ++j)
for(ii = 0; ii < nJ; ++ii)
if(J[ii][i]*J[ii][j])
Hg[i][j] = 1;

*nng = iData->nc * iData->nv;
return 0;
}

/*
* function calculates a jacobian matrix struct
Expand All @@ -375,7 +432,7 @@ static int local_jac_struct(IPOPT_DATA_ *iData, int * nng)
short **Hg;
int i,j,l,ii,nx, id;
int *cC,*lindex;
int nJ = (int)(iData->nx+iData->nc);
int nJ = (int)(iData->nJ);
id = 0;

J = iData->knowedJ;
Expand Down Expand Up @@ -559,12 +616,12 @@ static int optimizer_bounds_setings(DATA *data, IPOPT_DATA_ *iData)
{
int i, j, id;
for(i =0;i<iData->nx;++i){
check_nominal(iData, data->modelData.realVarsData[i].attribute.min, data->modelData.realVarsData[i].attribute.max, data->modelData.realVarsData[i].attribute.nominal, data->modelData.realVarsData[i].attribute.useNominal, i, fabs(iData->x0[i]));
iData->scalVar[i] = 1.0 / iData->vnom[i];
iData->scalf[i] = iData->scalVar[i];
check_nominal(iData, data->modelData.realVarsData[i].attribute.min, data->modelData.realVarsData[i].attribute.max, data->modelData.realVarsData[i].attribute.nominal, data->modelData.realVarsData[i].attribute.useNominal, i, fabs(iData->x0[i]));
iData->scalVar[i] = 1.0 / iData->vnom[i];
iData->scalf[i] = iData->scalVar[i];

iData->xmin[i] = data->modelData.realVarsData[i].attribute.min*iData->scalVar[i];
iData->xmax[i] = data->modelData.realVarsData[i].attribute.max*iData->scalVar[i];
iData->xmin[i] = data->modelData.realVarsData[i].attribute.min*iData->scalVar[i];
iData->xmax[i] = data->modelData.realVarsData[i].attribute.max*iData->scalVar[i];
}

iData->index_u = data->modelData.nVariablesReal - iData->nu;
Expand Down Expand Up @@ -643,7 +700,7 @@ static int optimizer_print_step(IPOPT_DATA_ *iData)
}


for(i=0 ,k=0,j =0; i<iData->NV; ++i,++j){
for(i = 0 ,k = 0,j = 0; i<iData->NV; ++i,++j){
if(j >= iData->nv){
j = 0; ++k;
}
Expand Down
11 changes: 0 additions & 11 deletions SimulationRuntime/c/optimization/mainOptimizer/ipoptODE.c
View file
Expand Up @@ -253,17 +253,6 @@ static int set_optimizer_flags(IPOPT_DATA_ *iData, IpoptProblem *nlp)
warningStreamPrint(LOG_STDOUT, 1, "not support ipopt_hesse=%s",cflags);
}

iData->useNumJac = 0;
cflags = (char*)omc_flagValue[FLAG_IPOPT_JAC];
if(cflags){
if(!strcmp(cflags,"NUM"))
iData->useNumJac = 1;
else if(!strcmp(cflags,"SYM") || !strcmp(cflags,"sym"))
iData->useNumJac = 0;
else
warningStreamPrint(LOG_STDOUT, 1, "not support ipopt_hesse=%s",cflags);
}

cflags = (char*)omc_flagValue[FLAG_LS_IPOPT];
if(cflags)
AddIpoptStrOption(*nlp, "linear_solver", cflags);
Expand Down
2 changes: 1 addition & 1 deletion SimulationRuntime/c/util/simulation_options.c
View file
Expand Up @@ -116,7 +116,7 @@ const char *FLAG_DETAILED_DESC[FLAG_MAX+1] = {
/* FLAG_INTERACTIVE */ "specify interactive simulation",
/* FLAG_IOM */ "value specifies the initialization optimization method",
/* FLAG_IPOPT_HESSE */ "value specifies the hessematrix for Ipopt(OMC, BFGS, const)",
/* FLAG_IPOPT_JAC */ "value specifies the jacobian for Ipopt(SYM, NUM)",
/* FLAG_IPOPT_JAC */ "value specifies the jacobian for Ipopt(SYM, NUM, NUMDENSE)",
/* FLAG_IPOPT_INIT */ "value specifies the initial guess for optimization (sim, const)",
/* FLAG_L */ "value specifies a time where the linearization of the model should be performed",
/* FLAG_LOG_FORMAT */ "value specifies the log format of the executable. -logFormat=text (default) or -logFormat=xml",
Expand Down

0 comments on commit 05c9498

Please sign in to comment.