Skip to content

Commit

Permalink
- fix under-determined initialization without initial equations
Browse files Browse the repository at this point in the history 
- fix error output for nonlinear systems

git-svn-id: https://openmodelica.org/svn/OpenModelica/trunk@10963 f25d12d1-65f4-0310-ae8a-bbce733d8d8e
  • Loading branch information
@lochel
lochel committed Jan 25, 2012
1 parent 71cfa4f commit d70a82f
Show file tree
Hide file tree
Showing 5 changed files with 43 additions and 49 deletions.
20 changes: 7 additions & 13 deletions Compiler/susan_codegen/SimCode/CodegenC.tpl
View file
Expand Up @@ -1934,7 +1934,7 @@ case SES_NONLINEAR(__) then
nls_xold[<%i0%>] = _<%namestr%>(1) /*old*/;
>>
;separator="\n"%>
solve_nonlinear_system(residualFunc<%index%>, -1, (void*)data);
solve_nonlinear_system(residualFunc<%index%>, data->modelData.equationInfo[SIM_PROF_EQ_<%index%>], (void*)data);
<%crefs |> name hasindex i0 => '<%cref(name)%> = nls_x[<%i0%>];' ;separator="\n"%>
end_nonlinear_system();<%inlineCrefs(context,crefs)%>
#ifdef _OMC_MEASURE_TIME
Expand Down Expand Up @@ -6646,7 +6646,7 @@ end literalExpConstArrayVal;
template equationInfo(list<SimEqSystem> eqs)
::=
match eqs
case {} then "const struct EQUATION_INFO equation_info[1] = {{0,NULL}};"
case {} then "const struct EQUATION_INFO equation_info[1] = {{0, NULL}};"
else
let &preBuf = buffer ""
let res =
Expand All @@ -6664,14 +6664,14 @@ template equationInfo(list<SimEqSystem> eqs)
//let var = '<%cref(componentRef)%>__varInfo'
//let &preBuf += 'const struct VAR_INFO *equationInfo_cref<%i0%> = &<%var%>;'
'"SES_ARRAY_CALL_ASSIGN <%i0%>",0,NULL'
case SES_ALGORITHM(__) then '"SES_ALGORITHM <%i0%>",0,NULL'
case SES_WHEN(__) then '"SES_WHEN <%i0%>",0,NULL'
case SES_LINEAR(__) then '"LINEAR<%index%>",0,NULL'
case SES_ALGORITHM(__) then '"SES_ALGORITHM <%i0%>", 0, NULL'
case SES_WHEN(__) then '"SES_WHEN <%i0%>", 0, NULL'
case SES_LINEAR(__) then '"LINEAR<%index%>", 0, NULL'
case SES_NONLINEAR(__) then
let &preBuf += 'const VAR_INFO** residualFunc<%index%>_crefs = (const VAR_INFO**)malloc(<%listLength(crefs)%>*sizeof(VAR_INFO*));<%\n%>'
let &preBuf += '<%crefs|>cr hasindex i0 => 'residualFunc<%index%>_crefs[<%i0%>] = &<%cref(cr)%>__varInfo;'; separator="\n"%>;'
'"residualFunc<%index%>",<%listLength(crefs)%>, residualFunc<%index%>_crefs'
case SES_MIXED(__) then '"MIXED<%index%>",0,NULL'
'"residualFunc<%index%>", <%listLength(crefs)%>, residualFunc<%index%>_crefs'
case SES_MIXED(__) then '"MIXED<%index%>", 0, NULL'
else '"unknown equation <%i0%>",0,NULL'%>}
>> ; separator=",\n"%>
};
Expand All @@ -6698,12 +6698,6 @@ template equationInfo(list<SimEqSystem> eqs)
; empty
%>
};
<% eqs |> eq hasindex i0 => match eq
case SES_MIXED(__)
case SES_LINEAR(__)
case SES_NONLINEAR(__) then '#define SIM_PROF_EQ_<%index%> <%i0%>'
; separator="\n"
%>
>>
end equationInfo;

Expand Down
36 changes: 18 additions & 18 deletions SimulationRuntime/c/math-support/initialization.c
View file
Expand Up @@ -159,15 +159,15 @@ static void computeInitialResidualScalingCoefficients(DATA *data, double nz, dou
DEBUG_INFO1(LOG_INIT, "initial residuals scaling coefficients [lambda=%g]", lambda);
for(j=0; j<data->modelData.nResiduals; ++j)
{
if(initialResidualScalingCoefficients[j] < 1e-42)
{
initialResidualScalingCoefficients[j] = 1.0;
DEBUG_INFO_AL2(LOG_INIT, " initial residual no. %d: %g [ineffective]", j, initialResidualScalingCoefficients[j]);
}
else
{
DEBUG_INFO_AL2(LOG_INIT, " initial residual no. %d: %g", j, initialResidualScalingCoefficients[j]);
}
if(initialResidualScalingCoefficients[j] < 1e-42)
{
initialResidualScalingCoefficients[j] = 1.0;
DEBUG_INFO_AL2(LOG_INIT, " initial residual no. %d: %g [ineffective]", j, initialResidualScalingCoefficients[j]);
}
else
{
DEBUG_INFO_AL2(LOG_INIT, " initial residual no. %d: %g", j, initialResidualScalingCoefficients[j]);
}
}

free(tmpInitialResidual1);
Expand Down Expand Up @@ -498,7 +498,7 @@ static int nelderMeadEx_initialization(DATA *data, long nz, double *z, char** zN
{
DEBUG_INFO1(LOG_INIT, "initialization-nr. %d", (int)l);

NelderMeadOptimization(nz, z, zNominal, initialResidualScalingCoefficients, lambda_step, STOPCR, NLOOP, DEBUG_FLAG(LOG_INIT) ? 100 : 0, &lambda, &iteration, leastSquareWithLambda, data, initialResiduals);
NelderMeadOptimization(nz, z, zNominal, initialResidualScalingCoefficients, lambda_step, STOPCR, NLOOP, DEBUG_FLAG(LOG_INIT) ? 10000 : 0, &lambda, &iteration, leastSquareWithLambda, data, initialResiduals);

storePreValues(data); /* save pre-values */
overwriteOldSimulationData(data); /* if there are non-linear equations */
Expand Down Expand Up @@ -619,10 +619,7 @@ static INIT_DATA *initializeInitData(DATA *data)
{
initData->zNominal[iz] = fabs(data->modelData.realVarsData[i].attribute.nominal);
if(initData->zNominal[iz] == 0.0)
{
DEBUG_INFO1(LOG_INIT, "set nominal(%s) := 1.0", data->modelData.realVarsData[i].info.name);
initData->zNominal[iz] = 1.0;
}
initData->z[iz] = data->modelData.realVarsData[i].attribute.start;
initData->zName[iz] = data->modelData.realVarsData[i].info.name;
iz++;
Expand All @@ -635,10 +632,7 @@ static INIT_DATA *initializeInitData(DATA *data)
{
initData->zNominal[iz] = fabs(data->modelData.realParameterData[i].attribute.nominal);
if(initData->zNominal[iz] == 0.0)
{
DEBUG_INFO1(LOG_INIT, "set nominal(%s) := 1.0", data->modelData.realParameterData[i].info.name);
initData->zNominal[iz] = 1.0;
}
initData->z[iz] = data->modelData.realParameterData[i].attribute.start;
initData->zName[iz] = data->modelData.realParameterData[i].info.name;
iz++;
Expand Down Expand Up @@ -669,7 +663,13 @@ static int initialize(DATA *data, int optiMethod)
/* no initial values to calculate. */
if(initData == NULL)
{
DEBUG_INFO(LOG_INIT, "no initial values to calculate");
DEBUG_INFO(LOG_INIT, "no variables to initialize");
return 0;
}

if(data->modelData.nInitEquations == 0)
{
DEBUG_INFO(LOG_INIT, "no initial equations");
return 0;
}

Expand Down Expand Up @@ -733,7 +733,7 @@ static int initialize(DATA *data, int optiMethod)

freeInitData(initData);
initData = initializeInitData(data);
/* no initial values to calculate. */
/* no initial values to calculate. (not possible to be here)*/
if(initData == NULL)
{
DEBUG_INFO(LOG_INIT, "no initial values to calculate");
Expand Down
22 changes: 11 additions & 11 deletions SimulationRuntime/c/math-support/matrix.h
View file
Expand Up @@ -85,7 +85,7 @@ void * _omc_hybrj_(void(*) (int*, double*, double*, double *, int*, int*, void*
printf("}\n"); \
} while(0)

#define solve_nonlinear_system_mixed(residual,no, userdata) do { \
#define solve_nonlinear_system_mixed(residual, no, userdata) do { \
int giveUp = 0; \
int retries = 0; \
int retries2 = 0; \
Expand Down Expand Up @@ -113,7 +113,7 @@ void * _omc_hybrj_(void(*) (int*, double*, double*, double *, int*, int*, void*
} else if (info >= 2 && info <= 5) { \
int i = 0; \
modelErrorCode=ERROR_NONLINSYS; \
DEBUG_INFO2(LOG_NONLIN_SYS,"error solving nonlinear system nr. %d at time %f",-1,data->localData[0]->timeValue); \
DEBUG_INFO2(LOG_NONLIN_SYS, "error solving nonlinear system nr. %d at time %f", no, data->localData[0]->timeValue); \
if (DEBUG_FLAG(LOG_NONLIN_SYS)) { \
for (i = 0; i < n; i++) { \
DEBUG_INFO_AL2(LOG_NONLIN_SYS," residual[%d] = %f",i,nls_fvec[i]); \
Expand All @@ -124,7 +124,7 @@ void * _omc_hybrj_(void(*) (int*, double*, double*, double *, int*, int*, void*
}\
} while(0) /* (no trailing ;)*/

#define solve_nonlinear_system(residual,no, userdata) do { \
#define solve_nonlinear_system(residual, no, userdata) do { \
int giveUp = 0; \
int retries = 0; \
int retries2 = 0; \
Expand All @@ -134,26 +134,26 @@ void * _omc_hybrj_(void(*) (int*, double*, double*, double *, int*, int*, void*
nls_diag,&mode,&factor,&nprint,&info,&nfev,nls_fjac,&ldfjac, \
nls_r,&lr,nls_qtf,nls_wa1,nls_wa2,nls_wa3,nls_wa4, userdata); \
if (info == 0) \
printErrorEqSyst(IMPROPER_INPUT,-1,data->localData[0]->timeValue); \
printErrorEqSyst(IMPROPER_INPUT, no, data->localData[0]->timeValue); \
if ((info == 4 || info == 5) && retries < 3) { /* first try to decrease factor*/ \
retries++; giveUp = 0; \
factor = factor / 10.0; \
if (DEBUG_FLAG(LOG_NONLIN_SYS)) \
printErrorEqSyst(NO_PROGRESS_FACTOR,-1,factor); \
printErrorEqSyst(NO_PROGRESS_FACTOR, no, factor); \
} else if ((info == 4 || info == 5) && retries < 5) { /* Then, try with different starting point*/ \
int i = 0; \
for (i = 0; i < n; i++) { nls_x[i] += 0.1; }; \
retries++; giveUp=0; \
if (DEBUG_FLAG(LOG_NONLIN_SYS)) \
printErrorEqSyst(NO_PROGRESS_START_POINT,no,1e-6); \
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++; giveUp = 0; \
int i = 0; \
for (i = 0; i < n; i++) { nls_x[i] = nls_xold[i]; } \
} else if (info >= 2 && info <= 5) { \
int i = 0; \
modelErrorCode=ERROR_NONLINSYS; \
printErrorEqSyst(ERROR_AT_TIME,-1,data->localData[0]->timeValue); \
printErrorEqSyst(ERROR_AT_TIME, no, data->localData[0]->timeValue); \
if (DEBUG_FLAG(LOG_NONLIN_SYS)) { \
for (i = 0; i < n; i++) { \
DEBUG_INFO_AL2(LOG_NONLIN_SYS," residual[%d] = %f",i,nls_fvec[i]); \
Expand All @@ -173,23 +173,23 @@ void * _omc_hybrj_(void(*) (int*, double*, double*, double *, int*, int*, void*
nls_diag,&mode,&factor,&nprint,&info,&nfev,&njev, \
nls_r,&lr,nls_qtf,nls_wa1,nls_wa2,nls_wa3,nls_wa4, userdata); \
if (info == 0) { \
printErrorEqSyst(IMPROPER_INPUT,-1,data->localData[0]->timeValue); \
printErrorEqSyst(IMPROPER_INPUT, no, data->localData[0]->timeValue); \
} \
if ((info == 4 || info == 5) && retries < 3) { /* First try to decrease factor*/ \
retries++; giveUp = 0; \
factor = factor / 10.0; \
if (sim_verbose) \
printErrorEqSyst(NO_PROGRESS_FACTOR,-1,factor); \
printErrorEqSyst(NO_PROGRESS_FACTOR, no, factor); \
} else if ((info == 4 || info == 5) && retries < 5) { /* Secondly, try with different starting point*/ \
int i = 0; \
for (i = 0; i < n; i++) { nls_x[i] += 0.1; }; \
retries++; giveUp=0; \
if (sim_verbose) \
printErrorEqSyst(NO_PROGRESS_START_POINT,-1,1e-6); \
printErrorEqSyst(NO_PROGRESS_START_POINT, no, 1e-6); \
} \
else if (info >= 2 && info <= 5) { \
modelErrorCode=ERROR_NONLINSYS; \
printErrorEqSyst(-1,data->localData[0]->timeValue); \
printErrorEqSyst(no, data->localData[0]->timeValue); \
} \
}\
} while(0) /* (no trailing ;)*/
Expand Down
12 changes: 6 additions & 6 deletions SimulationRuntime/c/util/varinfo.c
View file
Expand Up @@ -34,23 +34,23 @@
#include "varinfo.h"


void printErrorEqSyst(equationSystemError err, struct EQUATION_INFO eq, double var)
void printErrorEqSyst(equationSystemError err, EQUATION_INFO eq, double time)
{
switch (err) {
case ERROR_AT_TIME:
WARNING2("Error solving nonlinear system %s at time %g", eq.name, var);
WARNING2("Error solving nonlinear system %s at time %g", eq.name, time);
break;
case NO_PROGRESS_START_POINT:
WARNING2("Solving nonlinear system %s: iteration not making progress, trying with different starting points (+%g)", eq.name, var);
WARNING2("Solving nonlinear system %s: iteration not making progress, trying with different starting points (+%g)", eq.name, time);
break;
case NO_PROGRESS_FACTOR:
WARNING2("Solving nonlinear system %s: iteration not making progress, trying to decrease factor to %g", eq.name, var);
WARNING2("Solving nonlinear system %s: iteration not making progress, trying to decrease factor to %g", eq.name, time);
break;
case IMPROPER_INPUT:
WARNING2("improper input parameters to nonlinear eq. syst: %s at time %g", eq.name, var);
WARNING2("improper input parameters to nonlinear eq. syst: %s at time %g", eq.name, time);
break;
default:
WARNING3("Unknown equation system error: %d %s %g", err, eq.name, var);
WARNING3("Unknown equation system error: %d %s %g", err, eq.name, time);
break;
}
}
Expand Down
2 changes: 1 addition & 1 deletion SimulationRuntime/c/util/varinfo.h
View file
Expand Up @@ -38,7 +38,7 @@
extern "C" {
#endif

void printErrorEqSyst(equationSystemError, struct EQUATION_INFO, double var);
void printErrorEqSyst(equationSystemError err, EQUATION_INFO eq, double time);

void printInfo(FILE *stream, FILE_INFO info);

Expand Down

0 comments on commit d70a82f

Please sign in to comment.