Skip to content

Commit

Permalink
fix for parameters with nominal-attribute
Browse files Browse the repository at this point in the history 
git-svn-id: https://openmodelica.org/svn/OpenModelica/trunk@10641 f25d12d1-65f4-0310-ae8a-bbce733d8d8e
  • Loading branch information
@lochel
lochel committed Dec 5, 2011
1 parent cbbf705 commit 24f9b90
Show file tree
Hide file tree
Showing 3 changed files with 89 additions and 78 deletions.
20 changes: 12 additions & 8 deletions SimulationRuntime/c/math-support/events.c
View file
Expand Up @@ -260,21 +260,25 @@ void storeStartValues(_X_DATA* data)
SIMULATION_DATA *sData = data->localData[0];
MODEL_DATA *mData = &(data->modelData);

for(i=0; i<mData->nVariablesReal; ++i){
for(i=0; i<mData->nVariablesReal; ++i)
{
sData->realVars[i] = mData->realVarsData[i].attribute.start;
DEBUG_INFO2(LOG_DEBUG,"Set Real var %s = %g",mData->realVarsData[i].info.name, sData->realVars[i]);
DEBUG_INFO2(LOG_DEBUG, "Set Real var %s = %g", mData->realVarsData[i].info.name, sData->realVars[i]);
}
for(i=0; i<mData->nVariablesInteger; ++i){
for(i=0; i<mData->nVariablesInteger; ++i)
{
sData->integerVars[i] = mData->integerVarsData[i].attribute.start;
DEBUG_INFO2(LOG_DEBUG,"Set Integer var %s = %d",mData->integerVarsData[i].info.name, sData->integerVars[i]);
DEBUG_INFO2(LOG_DEBUG, "Set Integer var %s = %d", mData->integerVarsData[i].info.name, sData->integerVars[i]);
}
for(i=0; i<mData->nVariablesBoolean; ++i){
for(i=0; i<mData->nVariablesBoolean; ++i)
{
sData->booleanVars[i] = mData->booleanVarsData[i].attribute.start;
DEBUG_INFO2(LOG_DEBUG,"Set Boolean var %s = %s",mData->booleanVarsData[i].info.name, sData->booleanVars[i]?"true":"false");
DEBUG_INFO2(LOG_DEBUG, "Set Boolean var %s = %s", mData->booleanVarsData[i].info.name, sData->booleanVars[i]?"true":"false");
}
for(i=0; i<mData->nVariablesString; ++i){
for(i=0; i<mData->nVariablesString; ++i)
{
sData->stringVars[i] = copy_modelica_string((modelica_string_const)mData->stringVarsData[i].attribute.start);
DEBUG_INFO2(LOG_DEBUG,"Set String var %s = %s",mData->stringVarsData[i].info.name, sData->stringVars[i]);
DEBUG_INFO2(LOG_DEBUG, "Set String var %s = %s", mData->stringVarsData[i].info.name, sData->stringVars[i]);
}
}

Expand Down
99 changes: 51 additions & 48 deletions SimulationRuntime/c/math-support/initialization.c
View file
Expand Up @@ -105,28 +105,28 @@ const char *optiMethodStr[4] = {"unknown", "simplex", "nelder_mead_ex", "newuoa"
* \param z [in] vector of scaling-factors or NULL
* \param lambda [in]
*/
static double leastSquareWithLambda(long nz, double* z, double* scale, double lambda, _X_DATA* data, double* initialResiduals)
static double leastSquareWithLambda(long nz, double* z, double* zNominal, double lambda, _X_DATA* data, double* initialResiduals)
{
int indz = 0;
fortran_integer i = 0;
long j = 0;
double funcValue = 0;

for(i=0; i<data->modelData.nStates; i++)
for(i=0; i<data->modelData.nStates; ++i)
{
if(data->modelData.realVarsData[i].attribute.fixed==0)
{
data->localData[0]->realVars[i] = z[indz] * (scale ? scale[indz] : 1.0);
data->localData[0]->realVars[i] = z[indz] * (zNominal ? zNominal[indz] : 1.0);
indz++;
}
}

/* for real parameters */
for(i=0; i<data->modelData.nParametersReal; i++)
for(i=0; i<data->modelData.nParametersReal; ++i)
{
if(data->modelData.realParameterData[i].attribute.fixed == 0)
{
data->simulationInfo.realParameter[i] = z[indz] * (scale ? scale[indz] : 1.0);
data->simulationInfo.realParameter[i] = z[indz] * (zNominal ? zNominal[indz] : 1.0);
indz++;
}
}
Expand All @@ -136,8 +136,7 @@ static double leastSquareWithLambda(long nz, double* z, double* scale, double la
functionAlgebraics(data);

initial_residual(data, lambda, initialResiduals);

for(j=0; j<data->modelData.nResiduals; j++)
for(j=0; j<data->modelData.nResiduals; ++j)
{
funcValue += initialResiduals[j] * initialResiduals[j];
}
Expand All @@ -161,13 +160,13 @@ static void NelderMeadOptimization(long N,
double beta = 2; /* 1 < beta */
double gamma = 0.5; /* 0 < gamma < 1 */

double* simplex = (double*)calloc((N+1) * N,sizeof(double));
double* fvalues = (double*)calloc(N+1,sizeof(double));
double* simplex = (double*)calloc((N+1)*N, sizeof(double));
double* fvalues = (double*)calloc(N+1, sizeof(double));

double* xr = (double*)calloc(N,sizeof(double));
double* xe = (double*)calloc(N,sizeof(double));
double* xk = (double*)calloc(N,sizeof(double));
double* xbar = (double*)calloc(N,sizeof(double));
double* xr = (double*)calloc(N, sizeof(double));
double* xe = (double*)calloc(N, sizeof(double));
double* xk = (double*)calloc(N, sizeof(double));
double* xbar = (double*)calloc(N, sizeof(double));

double fxr;
double fxe;
Expand All @@ -180,7 +179,7 @@ static void NelderMeadOptimization(long N,
long x = 0;
long i = 0;

double lambda = 0.0; /* no lambda-control is activated */
double lambda = 0.0;
long iteration = 0;

/* check Memory */
Expand Down Expand Up @@ -394,36 +393,36 @@ static int reportResidualValue(double funcValue, _X_DATA* data, double* initialR
* nelderMead algorithm.
* This does not require a jacobian for the residuals.
*/
static int nelderMeadEx_initialization(_X_DATA *data, long nz, double *z, double *scale, double* initialResiduals)
static int nelderMeadEx_initialization(_X_DATA *data, long nz, double *z, double *zNominal, double* initialResiduals)
{
double STOPCR = 1.e-16;
double lambda_step = 0.1;
long NLOOP = 10000 * nz;

double funcValue = leastSquareWithLambda(nz, z, scale, 1.0, data, initialResiduals);
double funcValue = leastSquareWithLambda(nz, z, zNominal, 1.0, data, initialResiduals);

double lambda = 0;
long iteration = 0;

long l=0,i=0;
long l=0, i=0;

for(l=0; l<100 && funcValue > STOPCR; l++)
{
DEBUG_INFO1(LOG_INIT, "nelderMeadEx_initialization | initialization-nr. %d", (int)l);

/* down-scale */
for(i=0; i<nz; i++)
z[i] /= scale[i];
z[i] /= zNominal[i];

NelderMeadOptimization(nz, z, scale, lambda_step, STOPCR, NLOOP, useVerboseOutput(LOG_INIT) ? 100000 : 0, &lambda, &iteration, leastSquareWithLambda, data, initialResiduals);
NelderMeadOptimization(nz, z, zNominal, lambda_step, STOPCR, NLOOP, DEBUG_FLAG(LOG_INIT) ? 100000 : 0, &lambda, &iteration, leastSquareWithLambda, data, initialResiduals);

/* up-scale */
for(i=0; i<nz; i++)
z[i] *= scale[i];
z[i] *= zNominal[i];

if(DEBUG_FLAG(LOG_INIT))
{
INFO3("nelderMeadEx_initialization | iteration=%d / lambda=%g / f=%g", (int) iteration, lambda, leastSquareWithLambda(nz, z, scale, lambda, data, initialResiduals));
INFO3("nelderMeadEx_initialization | iteration=%d / lambda=%g / f=%g", (int) iteration, lambda, leastSquareWithLambda(nz, z, zNominal, lambda, data, initialResiduals));
for(i=0; i<nz; i++)
{
INFO_AL2("nelderMeadEx_initialization | states | %d: %g", (int) i, z[i]);
Expand All @@ -441,7 +440,7 @@ static int nelderMeadEx_initialization(_X_DATA *data, long nz, double *z, double
storePreValues(data);
overwriteOldSimulationData(data);

funcValue = leastSquareWithLambda(nz, z, scale, 1.0, data, initialResiduals);
funcValue = leastSquareWithLambda(nz, z, zNominal, 1.0, data, initialResiduals);
}

DEBUG_INFO1(LOG_INIT, "nelderMeadEx_initialization | leastSquare=%g", funcValue);
Expand All @@ -466,7 +465,7 @@ static int initialize(_X_DATA *data, int optiMethod)
long nz = 0;
int retVal = 0;
double *z = NULL;
double *scale = NULL;
double *zNominal = NULL;
double *initialResiduals = NULL;

DEBUG_INFO1(LOG_INIT, "initialization by method: %s", optiMethodStr[optiMethod]);
Expand All @@ -475,30 +474,18 @@ static int initialize(_X_DATA *data, int optiMethod)
DEBUG_INFO(LOG_INIT, "fixed attribute for states:");
for(i=0; i<data->modelData.nStates; ++i)
{
DEBUG_INFO2(LOG_INIT, "state %s(fixed=%s)", data->modelData.realVarsData[i].info.name, (data->modelData.realVarsData[i].attribute.fixed ? "true" : "false"));
if(data->modelData.realVarsData[i].attribute.fixed == 0)
{
DEBUG_INFO1(LOG_INIT, "state %s(fixed=false)", data->modelData.realVarsData[i].info.name);
++nz;
}
else
{
DEBUG_INFO1(LOG_INIT, "state %s(fixed=true)", data->modelData.realVarsData[i].info.name);
}
}

/* plus unfixed real-parameters */
DEBUG_INFO(LOG_INIT, "fixed attribute for parameters:");
for(i=0; i<data->modelData.nParametersReal; ++i)
{
DEBUG_INFO2(LOG_INIT, "parameter %s(fixed=%s)", data->modelData.realParameterData[i].info.name, (data->modelData.realParameterData[i].attribute.fixed ? "true" : "false"));
if(data->modelData.realParameterData[i].attribute.fixed == 0)
{
DEBUG_INFO1(LOG_INIT, "parameter %s(fixed=false)", data->modelData.realParameterData[i].info.name);
++nz;
}
else
{
DEBUG_INFO1(LOG_INIT, "parameter %s(fixed=true)", data->modelData.realParameterData[i].info.name);
}
}

DEBUG_INFO1(LOG_INIT, "number of non-fixed variables: %d", (int)nz);
Expand All @@ -511,17 +498,25 @@ static int initialize(_X_DATA *data, int optiMethod)
}

z = (double*)calloc(nz, sizeof(double));
scale = (double*)calloc(nz, sizeof(double));
zNominal = (double*)calloc(nz, sizeof(double));
initialResiduals = (double*) calloc(data->modelData.nResiduals, sizeof(double));
ASSERT(z, "out of memory");
ASSERT(scale, "out of memory");
ASSERT(zNominal, "out of memory");
ASSERT(initialResiduals, "out of memory");

/* fill z with the non-fixed variables from x and p */
for(i=0, iz=0; i<data->modelData.nStates; ++i)
{
if(data->modelData.realVarsData[i].attribute.fixed == 0)
{
scale[iz] = data->modelData.realVarsData[i].attribute.useNominal ? fabs(data->modelData.realVarsData[i].attribute.nominal) : 1;
zNominal[iz] = fabs(data->modelData.realVarsData[i].attribute.nominal);

if(zNominal[iz] == 0.0)
{
WARNING1("nominal for real parameter is zero > nominal(%s) := 1.0", data->modelData.realParameterData[i].info.name);
zNominal[iz] = 1.0;
}

z[iz++] = data->modelData.realVarsData[i].attribute.start;
}
}
Expand All @@ -531,14 +526,21 @@ static int initialize(_X_DATA *data, int optiMethod)
{
if(data->modelData.realParameterData[i].attribute.fixed == 0)
{
scale[iz] = data->modelData.realParameterData[i].attribute.useNominal ? fabs(data->modelData.realParameterData[i].attribute.nominal) : 1;
zNominal[iz] = fabs(data->modelData.realParameterData[i].attribute.nominal);

if(zNominal[iz] == 0.0)
{
WARNING1("nominal for real parameter is zero > nominal(%s) := 1.0", data->modelData.realParameterData[i].info.name);
zNominal[iz] = 1.0;
}

z[iz++] = data->modelData.realParameterData[i].attribute.start;
}
}

if(optiMethod == IOM_NELDER_MEAD_EX)
{
retVal = nelderMeadEx_initialization(data, nz, z, scale, initialResiduals);
retVal = nelderMeadEx_initialization(data, nz, z, zNominal, initialResiduals);
}
/*
else if(optiMethod == IOM_SIMPLEX)
Expand All @@ -553,12 +555,13 @@ static int initialize(_X_DATA *data, int optiMethod)
else
{
WARNING1("unrecognized option -iom %s", optiMethodStr[optiMethod]);
WARNING_AL("current options are: nelder_mead_ex, simplex or newuoa");
WARNING_AL("current options are: nelder_mead_ex"); /* WARNING_AL("current options are: nelder_mead_ex, simplex or newuoa"); */
retVal= -1;
}

free(z);
free(scale);
free(zNominal);
free(initialResiduals);
return retVal;
}

Expand All @@ -579,18 +582,18 @@ static int state_initialization(_X_DATA *data, int optiMethod)
storePreValues(data); /* to provide all valid pre-values */
overwriteOldSimulationData(data);

/* Initialize all relations that are ZeroCrossings */
/* initialize all relations that are ZeroCrossings */
bound_parameters(data);
update_DAEsystem(data);

/* And restore start values and helpvars */
resetAllHelpVars(data); /* ??? */
/* and restore start values and helpvars */
resetAllHelpVars(data);
storePreValues(data);

/* start with the real initialization */
data->simulationInfo.initial = 1; /* to evaluate when-equations with initial()-conditions */

retVal = initialize(data, IOM_NELDER_MEAD_EX);
retVal = initialize(data, optiMethod);

storeInitialValuesParam(data);
storePreValues(data); /* save pre-values */
Expand Down
48 changes: 26 additions & 22 deletions SimulationRuntime/c/simulation/simulation_input_xml.cpp
View file
Expand Up @@ -376,20 +376,20 @@ void read_input_xml(int argc, char **argv,
DEBUG_INFO2(LOG_DEBUG," read for %s readonly %d from init file",modelData->realVarsData[i].info.name,modelData->realVarsData[i].info.info.readonly);

/* read var attribute */
read_value(mi.rSta[i]["start"],&(modelData->realVarsData[i].attribute.start));
DEBUG_INFO2(LOG_DEBUG," read start-value for %s = %f from init file",modelData->realVarsData[i].info.name,modelData->realVarsData[i].attribute.start);
read_value(mi.rSta[i]["fixed"],(signed char*)&(modelData->realVarsData[i].attribute.fixed));
DEBUG_INFO2(LOG_DEBUG," read fixed for %s = %s from init file",modelData->realVarsData[i].info.name,(modelData->realVarsData[i].attribute.fixed)?"true":"false");
read_value(mi.rSta[i]["start"], &(modelData->realVarsData[i].attribute.start));
DEBUG_INFO2(LOG_DEBUG," read start-value for %s = %f from init file", modelData->realVarsData[i].info.name,modelData->realVarsData[i].attribute.start);
read_value(mi.rSta[i]["fixed"], (signed char*)&(modelData->realVarsData[i].attribute.fixed));
DEBUG_INFO2(LOG_DEBUG," read fixed for %s = %s from init file", modelData->realVarsData[i].info.name, (modelData->realVarsData[i].attribute.fixed)?"true":"false");
read_value(mi.rSta[i]["nominal"], &(modelData->realVarsData[i].attribute.nominal));
DEBUG_INFO2(LOG_DEBUG," read nominal-value for %s = %f from init file", modelData->realVarsData[i].info.name,modelData->realVarsData[i].attribute.nominal);
DEBUG_INFO2(LOG_DEBUG," read nominal-value for %s = %f from init file", modelData->realVarsData[i].info.name, modelData->realVarsData[i].attribute.nominal);

/* create a mapping for Alias variable to get the correct index */
mapAlias[(modelData->realVarsData[i].info.name)]=i;
}

/* Read stateDerivatives static data */
for(int i = 0; i < modelData->nStates; i++) {

for(int i = 0; i < modelData->nStates; i++)
{
DEBUG_INFO(LOG_DEBUG, "Read xml file for stateDerivatives");
/* read var info */
read_value(mi.rDer[i]["name"], &(modelData->realVarsData[modelData->nStates + i].info.name));
Expand All @@ -412,18 +412,20 @@ void read_input_xml(int argc, char **argv,
DEBUG_INFO2(LOG_DEBUG," read for %s readonly %d from init file",modelData->realVarsData[modelData->nStates+i].info.name,modelData->realVarsData[modelData->nStates+i].info.info.readonly);

/* read var attribute */
read_value(mi.rDer[i]["start"],&(modelData->realVarsData[modelData->nStates+i].attribute.start));
DEBUG_INFO2(LOG_SOLVER," read start-value for %s = %f from init file",modelData->realVarsData[modelData->nStates+i].info.name,modelData->realVarsData[modelData->nStates+i].attribute.start);
read_value(mi.rDer[i]["fixed"],(signed char*)&(modelData->realVarsData[modelData->nStates+i].attribute.fixed));
DEBUG_INFO2(LOG_SOLVER," read fixed for %s = %s from init file",modelData->realVarsData[modelData->nStates+i].info.name,(modelData->realVarsData[modelData->nStates+i].attribute.fixed)?"true":"false");
read_value(mi.rDer[i]["start"], &(modelData->realVarsData[modelData->nStates+i].attribute.start));
DEBUG_INFO2(LOG_SOLVER," read start-value for %s = %f from init file", modelData->realVarsData[modelData->nStates+i].info.name, modelData->realVarsData[modelData->nStates+i].attribute.start);
read_value(mi.rDer[i]["fixed"], (signed char*)&(modelData->realVarsData[modelData->nStates+i].attribute.fixed));
DEBUG_INFO2(LOG_SOLVER," read fixed for %s = %s from init file", modelData->realVarsData[modelData->nStates+i].info.name, (modelData->realVarsData[modelData->nStates+i].attribute.fixed)?"true":"false");
read_value(mi.rSta[i]["nominal"], &(modelData->realVarsData[modelData->nStates+i].attribute.nominal));
DEBUG_INFO2(LOG_DEBUG," read nominal-value for %s = %f from init file", modelData->realVarsData[i].info.name, modelData->realVarsData[modelData->nStates+i].attribute.nominal);

/* create a mapping for Alias variable to get the correct index */
mapAlias[(modelData->realVarsData[modelData->nStates+i].info.name)]= modelData->nStates+i;
}

/* Read real algebraics static data */
for(int i = 0; i < (modelData->nVariablesReal - 2*modelData->nStates); i++) {

for(int i = 0; i < (modelData->nVariablesReal - 2*modelData->nStates); i++)
{
int j = 2*modelData->nStates + i;
DEBUG_INFO(LOG_DEBUG, "Read xml file for real algebracis");
/* read var info */
Expand Down Expand Up @@ -457,8 +459,8 @@ void read_input_xml(int argc, char **argv,
}

/* Read integer variables static data */
for(int i = 0; i < modelData->nVariablesInteger; i++) {

for(int i = 0; i < modelData->nVariablesInteger; i++)
{
DEBUG_INFO(LOG_DEBUG, "Read xml file for integer algebracis");
/* read var info */
read_value(mi.iAlg[i]["name"], &(modelData->integerVarsData[i].info.name));
Expand Down Expand Up @@ -491,8 +493,8 @@ void read_input_xml(int argc, char **argv,
}

/* Read boolean variables static data */
for(int i = 0; i < modelData->nVariablesBoolean; i++) {

for(int i = 0; i < modelData->nVariablesBoolean; i++)
{
DEBUG_INFO(LOG_DEBUG, "Read xml file for boolean algebracis");
/* read var info */
read_value(mi.bAlg[i]["name"], &(modelData->booleanVarsData[i].info.name));
Expand Down Expand Up @@ -560,8 +562,8 @@ void read_input_xml(int argc, char **argv,
* Real all parameter
*/
/* Read Parameters static data */
for(int i = 0; i < modelData->nParametersReal; i++) {

for(int i = 0; i < modelData->nParametersReal; i++)
{
DEBUG_INFO(LOG_DEBUG, "Read xml file for real parameters");
/* read var info */
read_value(mi.rPar[i]["name"], &(modelData->realParameterData[i].info.name));
Expand All @@ -584,10 +586,12 @@ void read_input_xml(int argc, char **argv,
DEBUG_INFO2(LOG_DEBUG," read for %s readonly %d from init file",modelData->realParameterData[i].info.name,modelData->realParameterData[i].info.info.readonly);

/* read var attribute */
read_value(mi.rPar[i]["start"],&(modelData->realParameterData[i].attribute.start));
DEBUG_INFO2(LOG_SOLVER," read start-value for %s = %f from init file",modelData->realParameterData[i].info.name,modelData->realParameterData[i].attribute.start);
read_value(mi.rPar[i]["start"], &(modelData->realParameterData[i].attribute.start));
DEBUG_INFO2(LOG_SOLVER, "read start-value for %s = %f from init file", modelData->realParameterData[i].info.name, modelData->realParameterData[i].attribute.start);
read_value(mi.rPar[i]["fixed"],(signed char*)&(modelData->realParameterData[i].attribute.fixed));
DEBUG_INFO2(LOG_SOLVER," read fixed for %s = %s from init file",modelData->realParameterData[i].info.name,modelData->realParameterData[i].attribute.fixed?"true":"false");
DEBUG_INFO2(LOG_SOLVER, "read fixed for %s = %s from init file", modelData->realParameterData[i].info.name, modelData->realParameterData[i].attribute.fixed?"true":"false");
read_value(mi.rPar[i]["nominal"], &(modelData->realParameterData[i].attribute.nominal));
DEBUG_INFO2(LOG_SOLVER, "read nominal for %s = %s from init file", modelData->realParameterData[i].info.name, modelData->realParameterData[i].attribute.nominal);

/* create a mapping for Alias variable to get the correct index */
mapAliasParam[(modelData->realParameterData[i].info.name)]=i;
Expand Down

0 comments on commit 24f9b90

Please sign in to comment.