Skip to content

Commit

Permalink
- added bugfix for parameter initialization and activated the testcase
Browse files Browse the repository at this point in the history
 - added and changed some debug flag in matrix.h for LOG_NONLINSYS


git-svn-id: https://openmodelica.org/svn/OpenModelica/trunk@8506 f25d12d1-65f4-0310-ae8a-bbce733d8d8e
  • Loading branch information
Willi Braun committed Apr 5, 2011
1 parent 4db34f5 commit 202d14e
Show file tree
Hide file tree
Showing 3 changed files with 89 additions and 101 deletions.
2 changes: 1 addition & 1 deletion Compiler/susan_codegen/SimCode/SimCodeC.tpl
Expand Up @@ -290,7 +290,7 @@ case MODELINFO(varInfo=VARINFO(__), vars=SIMVARS(__)) then
memcpy(data->stringVariables.alias,omc__stringAlias,sizeof(DATA_STRING_ALIAS)*data->stringVariables.nAlias);
};

static char init_fixed[NX+NX+NY+NYINT+NYBOOL+NYSTR+NP+NPINT+NPBOOL+NPSTR] = {
static char init_fixed[NX+NX+NY+NYINT+NYBOOL+NP+NPINT+NPBOOL] = {
<%{(vars.stateVars |> SIMVAR(__) =>
'<%globalDataFixedInt(isFixed)%> /* <%crefStr(name)%> */'
;separator=",\n"),
Expand Down
153 changes: 82 additions & 71 deletions c_runtime/matrix.h
Expand Up @@ -96,36 +96,39 @@ void * hybrj_(void(*) (int *,double*,double*,double *,int*, int*),
nls_diag,&mode,&factor,&nprint,&info,&nfev,nls_fjac,&ldfjac, \
nls_r,&lr,nls_qtf,nls_wa1,nls_wa2,nls_wa3,nls_wa4); \
if (info == 0) { \
printf("improper input parameters to nonlinear eq. syst %s:%d.\n", __FILE__, __LINE__); \
if (sim_verbose >= LOG_NONLIN_SYS) { \
printf("improper input parameters to nonlinear eq. syst %s:%d.\n", __FILE__, __LINE__); \
} \
} \
if ((info == 4 || info == 5 )&& retries < 3) { /* first try to decrease factor*/ \
retries++; giveUp = 0; \
factor = factor / 10.0; \
if (sim_verbose & LOG_NONLIN_SYS) \
if ((info == 4 || info == 5 )&& retries < 3) { /* first try to decrease factor*/ \
retries++; giveUp = 0; \
factor = factor / 10.0; \
if (sim_verbose >= LOG_NONLIN_SYS) { \
printf("Solving nonlinear system: iteration not making progress, trying to decrease factor to %f\n",factor); \
} else if ((info == 4 || info == 5) && retries < 5) { /* Then, try with different starting point*/ \
int i; \
for (i=0; i < n; i++) { nls_x[i]+=0.1; }; \
retries++; giveUp=0; \
if (sim_verbose & LOG_NONLIN_SYS) \
retries++; giveUp=0; \
if (sim_verbose >= LOG_NONLIN_SYS) { \
printf("Solving nonlinear system: iteration not making progress, trying with different starting points (+1e-6)\n"); \
} else if ((info == 4 || info == 5) && retries2 < 1) { /*Then try with old values (instead of extrapolating )*/ \
retries = 0; retries2++; \
int i; \
for (i=0; i < n; i++) { nls_x[i] = nls_xold[i]; } \
} \
else if (info >= 2 && info <= 5) { \
int i; \
modelErrorCode=ERROR_NONLINSYS; \
printf("error solving nonlinear system nr. %d at time %f\n",no,time); \
if (sim_verbose & LOG_NONLIN_SYS) { \
for (i = 0; i<n; i++) { \
printf(" residual[%d] = %f\n",i,nls_fvec[i]); \
printf(" x[%d] = %f\n",i,nls_x[i]); \
} \
} else if ((info == 4 || info == 5) && retries2 < 1) { /*Then try with old values (instead of extrapolating )*/ \
retries = 0; retries2++; \
int i; \
for (i=0; i < n; i++) { nls_x[i] = nls_xold[i]; } \
} else if (info >= 2 && info <= 5) { \
int i; \
modelErrorCode=ERROR_NONLINSYS; \
if (sim_verbose >= LOG_NONLIN_SYS) { \
printf("error solving nonlinear system nr. %d at time %f\n",no,time); \
} \
if (sim_verbose >= LOG_NONLIN_SYS) { \
for (i = 0; i<n; i++) { \
printf(" residual[%d] = %f\n",i,nls_fvec[i]); \
printf(" x[%d] = %f\n",i,nls_x[i]); \
} \
} \
} \
} \
}\
}\
} while(0) /* (no trailing ;)*/

#define solve_nonlinear_system(residual,no) do { \
Expand All @@ -149,18 +152,18 @@ void * hybrj_(void(*) (int *,double*,double*,double *,int*, int*),
int i; \
for (i=0; i < n; i++) { nls_x[i]+=0.1; }; \
retries++; giveUp=0; \
if (sim_verbose & LOG_NONLIN_SYS) \
if (sim_verbose >= LOG_NONLIN_SYS) \
printErrorEqSyst(NO_PROGRESS_START_POINT,no,1e-6); \
} else if ((info == 4 || info == 5) && retries2 < 1) { /*Then try with old values (instead of extrapolating )*/ \
retries = 0; retries2++; \
} else if ((info == 4 || info == 5) && retries2 < 1) { /*Then try with old values (instead of extrapolating )*/ \
retries = 0; retries2++; \
int i; \
for (i=0; i < n; i++) { nls_x[i] = nls_xold[i]; } \
} \
else if (info >= 2 && info <= 5) { \
int i; \
modelErrorCode=ERROR_NONLINSYS; \
printErrorEqSyst(ERROR_AT_TIME,no,time); \
if (sim_verbose & LOG_NONLIN_SYS) { \
if (sim_verbose >= LOG_NONLIN_SYS) { \
for (i = 0; i<n; i++) { \
printf(" residual[%d] = %f\n",i,nls_fvec[i]); \
printf(" x[%d] = %f\n",i,nls_x[i]); \
Expand Down Expand Up @@ -192,12 +195,12 @@ void * hybrj_(void(*) (int *,double*,double*,double *,int*, int*),
retries++; giveUp=0; \
if (sim_verbose) \
printErrorEqSyst(NO_PROGRESS_START_POINT,no,1e-6); \
} \
else if (info >= 2 && info <= 5) { \
modelErrorCode=ERROR_NONLINSYS; \
printErrorEqSyst(no,time); \
} \
}\
} \
else if (info >= 2 && info <= 5) { \
modelErrorCode=ERROR_NONLINSYS; \
printErrorEqSyst(no,time); \
} \
}\
} while(0) /* (no trailing ;)*/

#define declare_matrix(A,nrows,ncols) double *A = real_alloc(nrows*ncols); \
Expand Down Expand Up @@ -225,10 +228,12 @@ for(int i=0; i<n; i++) ipiv[i] = 0; \
integer info; /* output */ \
dgesv_(&n,&nrhs,&A[0],&lda,ipiv,&b[0],&ldb,&info); \
if (info < 0) { \
printf("Error solving linear system of equations (no. %d) at time %f. Argument %d illegal.\n",id,localData->timeValue,info); \
if (sim_verbose >= LOG_NONLIN_SYS) \
printf("Error solving linear system of equations (no. %d) at time %f. Argument %d illegal.\n",id,localData->timeValue,info); \
} \
else if (info > 0) { \
printf("Error solving linear system of equations (no. %d) at time %f, system is singular.\n",id,localData->timeValue); \
if (sim_verbose >= LOG_NONLIN_SYS) \
printf("Error solving linear system of equations (no. %d) at time %f, system is singular.\n",id,localData->timeValue); \
} \
delete [] ipiv; \
} while (0) /* (no trailing ; ) */
Expand All @@ -242,7 +247,8 @@ for(int i=0; i<n; i++) ipiv[i] = 0; \
integer info; /* output */ \
dgesv_(&n,&nrhs,&A[0],&lda,ipiv,&b[0],&ldb,&info); \
if (info < 0) { \
printf("Error solving linear system of equations (no. %d) at time %f. Argument %d illegal.\n",id,localData->timeValue,info); \
if (sim_verbose >= LOG_NONLIN_SYS) \
printf("Error solving linear system of equations (no. %d) at time %f. Argument %d illegal.\n",id,localData->timeValue,info); \
} \
else if (info > 0) { \
found_solution=-1; \
Expand Down Expand Up @@ -302,47 +308,52 @@ int ldfjac = size;
((($P$old##v)-($P$old2##v))/(localData->oldTime-localData->oldTime2)*localData->timeValue \
+(localData->oldTime*($P$old2##v)-localData->oldTime2*($P$old##v))/ \
(localData->oldTime-localData->oldTime2))

#define mixed_equation_system(size) do { \
int found_solution = 0; \
int cur_value_indx=0; \
do { \
double discrete_loc[size]; \
double discrete_loc2[size]
int found_solution = 0; \
int cur_value_indx=0; \
do { \
double discrete_loc[size]; \
double discrete_loc2[size]

#define mixed_equation_system_end(size) } while (!found_solution); \
} while(0)

#define check_discrete_values(size,numValues) do {int i; \
if (found_solution == -1) { /*system of equations failed*/ \
found_solution=0; \
} else { \
found_solution = 1; \
for (i=0; i < size; i++) { \
if (fabs((discrete_loc[i] - discrete_loc2[i])) > 1e-12) {\
found_solution=0;\
}\
}\
}\
if (!found_solution ) { \
cur_value_indx++; \
if (cur_value_indx >= numValues/size) { \
found_solution=-1; \
} else {\
/* try next set of values*/ \
for (i=0; i < size; i++) { \
*loc_ptrs[i]=values[cur_value_indx*size+i]; \
} \
} \
} \
if (found_solution && sim_verbose) { /* we found a solution*/ \
{int i; \
printf("Result of mixed system discrete variables:\n"); \
for (i=0;i<size;i++) { \
printf("%s = %f pre(%s)= %f\n",getNameReal(loc_ptrs[i]),*loc_ptrs[i], \
getNameReal(loc_ptrs[i]),pre(*loc_ptrs[i])); \
} \
#define check_discrete_values(size,numValues) \
do { \
int i; \
if (found_solution == -1) { \
/*system of equations failed */ \
found_solution=0; \
} else { \
found_solution = 1; \
for (i=0; i < size; i++) { \
if (fabs((discrete_loc[i] - discrete_loc2[i])) > 1e-12) {\
found_solution=0;\
}\
}\
}\
if (!found_solution ) { \
cur_value_indx++; \
if (cur_value_indx >= numValues/size) { \
found_solution=-1; \
} else {\
/* try next set of values*/ \
for (i=0; i < size; i++) { \
*loc_ptrs[i]=values[cur_value_indx*size+i]; \
} \
} \
} \
/* we found a solution*/ \
if (found_solution && (sim_verbose >= LOG_NONLIN_SYS)){ \
int i; \
printf("Result of mixed system discrete variables:\n"); \
for (i=0;i<size;i++) { \
printf("%s = %f pre(%s)= %f\n",getNameReal(loc_ptrs[i]),*loc_ptrs[i], \
getNameReal(loc_ptrs[i]),pre(*loc_ptrs[i])); \
} \
} \
} \
} while(0)


#endif
35 changes: 6 additions & 29 deletions c_runtime/simulation_init.cpp
Expand Up @@ -40,26 +40,17 @@ void leastSquare(long *nz, double *z, double *funcValue)
{
int ind, indAct, indz;
int startIndPar = 2*globalData->nStates+globalData->nAlgebraic+globalData->intVariables.nAlgebraic+globalData->boolVariables.nAlgebraic;
for (ind=0, indAct=0, indz=0; ind<globalData->nStates; ind++)

for (ind=0, indAct=0, indz=0; ind<globalData->nStates; ind++)
if (globalData->initFixed[indAct++]==0)
globalData->states[ind] = z[indz++];

// for real parameters
for (ind=0,indAct=startIndPar; ind<globalData->nParameters; ind++)
if (globalData->initFixed[indAct++]==0)
globalData->parameters[ind] = z[indz++];

// for int parameters
for (ind=0,indAct=startIndPar+globalData->nParameters; ind<globalData->intVariables.nParameters; ind++) {
if (globalData->initFixed[indAct++]==0)
globalData->intVariables.parameters[ind] = (modelica_integer)z[indz++];
}

// for bool parameters
for (ind=0,indAct=startIndPar+globalData->nParameters+globalData->intVariables.nParameters; ind<globalData->boolVariables.nParameters; ind++) {
if (globalData->initFixed[indAct++]==0)
globalData->boolVariables.parameters[ind] = (modelica_boolean)z[indz++];
}

bound_parameters();
functionODE();
functionAlgebraics();

Expand Down Expand Up @@ -232,16 +223,14 @@ int initialize(const std::string*method)

nz++;
}

}

int startIndPar = 2*globalData->nStates+globalData->nAlgebraic+globalData->intVariables.nAlgebraic+globalData->boolVariables.nAlgebraic;
int endIndPar = startIndPar+globalData->nParameters+globalData->intVariables.nParameters+globalData->boolVariables.nParameters;
int endIndPar = startIndPar+globalData->nParameters;
for (ind=startIndPar;ind<endIndPar; ind++){
if (globalData->initFixed[ind]==0){
if (sim_verbose >= LOG_INIT)
printf("Variable %s is unfixed.\n",globalData->statesNames[ind].name);

printf("Variable %s is unfixed.\n",globalData->parametersNames[ind-startIndPar].name);
nz++;
}
}
Expand Down Expand Up @@ -277,18 +266,6 @@ int initialize(const std::string*method)
z[indz++] = globalData->parameters[ind];
}

// for int parameters
for (ind=0,indAct=startIndPar+globalData->nParameters; ind<globalData->intVariables.nParameters; ind++) {
if (globalData->initFixed[indAct++]==0)
z[indz++] = globalData->intVariables.parameters[ind];
}

// for bool parameters
for (ind=0,indAct=startIndPar+globalData->nParameters+globalData->intVariables.nParameters; ind<globalData->boolVariables.nParameters; ind++) {
if (globalData->initFixed[indAct++]==0)
z[indz++] = globalData->boolVariables.parameters[ind];
}

int retVal=0;
if (init_method == std::string("simplex")) {
retVal = simplex_initialization(nz,z);
Expand Down

0 comments on commit 202d14e

Please sign in to comment.