# Changeset 8512

Ignore:
Timestamp:
Mar 1, 2016 9:15:21 AM (4 years ago)
Message:

#3879 Added a check for a valid minimum norm solution in case of a singular Jacobian.
#4475 In case we have an entire row in the Jacobian equal to 0, use heuristic residual scaling instead of the maximum.

Location:
trunk
Files:
9 edited

Unmodified
Removed
• ## trunk/Python/src/tests_jmodelica/files/Modelica/Singular.mo

 r8065 (sin(x*y*z)+z) = 0; end NonLinear5; model NoMinimumNormSolution Real x,y,z,v; parameter Real a11 = 1; parameter Real a12 = 2; parameter Real a13 = 3; parameter Real a21 = 1; parameter Real a22 = 1; parameter Real a23 = 1; Real a31(start=0); Real a32(start=0); Real a33(start=0); parameter Real b[3] = {1,2,3}; equation a33=a31+a32; a33=a31^2+a32^2; a33=a31-a32; a11*x+a12*y+a13*z = b[1]; a21*x+a22*y+a23*z = b[2]; a31*x+a32*y+a33*z = b[3]; der(v) = time; end NoMinimumNormSolution; end Singular;
• ## trunk/Python/src/tests_jmodelica/simulation/test_assimulo_interface_fmi.py

 r8507 compile_fmu("Singular.NonLinear4", file_name) compile_fmu("Singular.NonLinear5", file_name) compile_fmu("Singular.NoMinimumNormSolution", file_name) @testattr(stddist = True) nose.tools.assert_almost_equal(res["z"][0] ,0.000000000) nose.tools.assert_almost_equal(res["z"][-1] ,-1.000000000) @testattr(stddist = True) def test_no_valid_minimum_norm_sol(self): model = load_fmu("NoMinimumNormSolution.fmu", log_level=3) model.set("_log_level", 3) model.set_log_level(log_level) nose.tools.assert_raises(FMUException, model.initialize) class Test_FMI_ODE_CS_2:
• ## trunk/RuntimeLibrary/src/jmi/jmi_block_solver.c

 r8507 #include #include #include #include #include #include #include #include #include "jmi_log.h" jmi_block_solver_options_t* options, void* problem_data){ int i; jmi_block_solver_t* block_solver = (jmi_block_solver_t*)calloc(1, sizeof(jmi_block_solver_t)); if(!block_solver) { block_solver->dx=(jmi_real_t*)calloc(n,sizeof(jmi_real_t));;                /**< \brief Work vector for the seed vector */ block_solver->f_scale = N_VNew_Serial(n); block_solver->scale_update_time = -1.0; if(n>0) { block_solver->J = NewDenseMat(n ,n); block_solver->J_scale = NewDenseMat(n ,n); } block_solver->res = (jmi_real_t*)calloc(n,sizeof(jmi_real_t)); block_solver->max = (jmi_real_t*)calloc(n,sizeof(jmi_real_t)); block_solver->nominal = (jmi_real_t*)calloc(n,sizeof(jmi_real_t)); block_solver->residual_nominal = (jmi_real_t*)calloc(n+1,sizeof(jmi_real_t)); block_solver->residual_heuristic_nominal = (jmi_real_t*)calloc(n+1,sizeof(jmi_real_t)); block_solver->initial = (jmi_real_t*)calloc(n,sizeof(jmi_real_t)); block_solver->start_set = (jmi_real_t*)calloc(n,sizeof(jmi_real_t)); block_solver->cur_time = 0; block_solver->using_max_min_scaling_flag = 0; /* Not using max/min scaling */ /* Initial equation scaling is 1.0 */ N_VConst_Serial(1.0,block_solver->f_scale); /* Make sure that residual_nominals and heuristic nominals always have values */ for(i=0; iresidual_nominal[i] = 1.0; block_solver->residual_heuristic_nominal[i] = 1.0; } switch(options->solver) { case JMI_KINSOL_SOLVER: { free(block_solver->dx); N_VDestroy_Serial(block_solver->f_scale); if(block_solver->n > 0) { DestroyMat(block_solver->J); DestroyMat(block_solver->J_scale); } free(block_solver->res); free(block_solver->dres); free(block_solver->max); free(block_solver->nominal); free(block_solver->residual_nominal); free(block_solver->residual_heuristic_nominal); free(block_solver->initial); free(block_solver->value_references); block_solver->F(block_solver->problem_data,real_vrs,block_solver->res,JMI_BLOCK_VALUE_REFERENCE); block_solver->F(block_solver->problem_data,block_solver->start_set,block_solver->res,JMI_BLOCK_START_SET); jmi_setup_f_residual_scaling(block_solver); for (i=0;in;i++) { } /* Compute appropriate equation scaling and function tolerance based on Jacobian J, nominal values (block->nominal) and current point (block->x). Store result in block->f_scale. */ void jmi_update_f_scale(jmi_block_solver_t *block) { int i, N = block->n; realtype curtime = block->cur_time; realtype* scale_ptr = 0; jmi_block_solver_options_t* bsop = block->options; int use_scaling_flag = bsop->residual_equation_scaling_mode; realtype tol = bsop->res_tol; block->scale_update_time = curtime; block->force_rescaling = 0; if(bsop->residual_equation_scaling_mode != jmi_residual_scaling_none) { /* Zero out the scales initially if we're modify this. */ N_VConst_Serial(0,block->f_scale); } scale_ptr = N_VGetArrayPointer(block->f_scale); /* Form scaled Jacobian as needed for automatic scaling and condition number checking*/ if (bsop->residual_equation_scaling_mode == jmi_residual_scaling_auto || bsop->residual_equation_scaling_mode == jmi_residual_scaling_aggressive_auto || bsop->residual_equation_scaling_mode == jmi_residual_scaling_full_jacobian_auto || bsop->residual_equation_scaling_mode == jmi_residual_scaling_hybrid || bsop->residual_equation_scaling_mode == jmi_residual_scaling_manual) { for (i = 0; i < N; i++) { int j; /* column scaling is formed by max(abs(nominal), abs(actual_value)) */ realtype xscale = RAbs(block->nominal[i]); realtype x = RAbs(block->x[i]); realtype* col_ptr; realtype* scaled_col_ptr; if(x < xscale) x = xscale; if(x < tol) x = tol; col_ptr = DENSE_COL(block->J, i); scaled_col_ptr = DENSE_COL(block->J_scale, i); /* row scaling is product of Jac entry and column scaling */ for (j = 0; j < N; j++) { realtype dres = col_ptr[j]; realtype fscale; fscale = dres * x; scaled_col_ptr[j] = fscale; scale_ptr[j] = MAX(scale_ptr[j], RAbs(fscale)); } } } /* put in manual non-zero scales in hybrid & manual cases */ if(    (bsop->residual_equation_scaling_mode == jmi_residual_scaling_manual) || (bsop->residual_equation_scaling_mode == jmi_residual_scaling_hybrid))    { for(i = 0; i < N; i++) { if(block->residual_nominal[i] != 0) { scale_ptr[i] = block->residual_nominal[i]; } } } if(bsop->residual_equation_scaling_mode == jmi_residual_scaling_auto || bsop->residual_equation_scaling_mode == jmi_residual_scaling_aggressive_auto || bsop->residual_equation_scaling_mode == jmi_residual_scaling_full_jacobian_auto || bsop->residual_equation_scaling_mode == jmi_residual_scaling_hybrid) { block->using_max_min_scaling_flag = 0; /* NOT using max/min scaling */ /* check that scaling factors have reasonable magnitude */ for(i = 0; i < N; i++) { if(scale_ptr[i] == 0.0) { jmi_log_node(block->log, logInfo, "JacobianZeroRow", "The derivatives for are all equal 0. Using heuristic residual nominals for scaling in ", i, block->label); scale_ptr[i] = block->residual_heuristic_nominal[i]; } else if(scale_ptr[i] < 1/bsop->max_residual_scaling_factor) { int j, maxInJacIndex = 0; realtype **jac = block->J_scale->cols; realtype maxInJac = RAbs(jac[0][i]); block->using_max_min_scaling_flag = 1; /* Using maximum scaling, singular Jacobian? */ for(j = 0; j < N; j++) { if(RAbs(maxInJac) < RAbs(jac[j][i])) { maxInJac = jac[j][i]; maxInJacIndex = j; } } scale_ptr[i] = bsop->max_residual_scaling_factor; jmi_log_node(block->log, logWarning, "MaxScalingUsed", "Poor scaling in . " "Consider changes in the model. Partial derivative of with respect to is ().", block->label, i, block->value_references[maxInJacIndex], DENSE_ELEM(block->J, i, maxInJacIndex) ,maxInJac); } else if(scale_ptr[i] > 1/bsop->min_residual_scaling_factor) { int j, maxInJacIndex = 0; realtype **jac = block->J_scale->cols; realtype maxInJac = RAbs(jac[0][i]); /* Likely not a problem: solver->using_max_min_scaling_flag = 1; -- Using minimum scaling */ for(j = 0; j < N; j++) { if(RAbs(maxInJac) < RAbs(jac[j][i])) { maxInJac = jac[j][i]; maxInJacIndex = j; } } if (scale_ptr[i] > block->residual_heuristic_nominal[i]/bsop->min_residual_scaling_factor) { scale_ptr[i] = bsop->min_residual_scaling_factor/block->residual_heuristic_nominal[i]; jmi_log_node(block->log, logWarning, "MinScalingUsed", "Poor scaling in . " "Consider changes in the model. Partial derivative of with respect to is ().", block->label, i, block->value_references[maxInJacIndex], DENSE_ELEM(block->J, i, maxInJacIndex) ,maxInJac); } else { scale_ptr[i] = 1/scale_ptr[i]; jmi_log_node(block->log, logWarning, "MinScalingExceeded", "Poor scaling in but will allow it based on structure of block. " "Consider changes in the model. Partial derivative of with respect to is ().", block->label, i, block->value_references[maxInJacIndex], DENSE_ELEM(block->J, i, maxInJacIndex) ,maxInJac); } } else scale_ptr[i] = 1/scale_ptr[i]; } if (block->callbacks->log_options.log_level >= 4) { jmi_log_node_t outer = jmi_log_enter_fmt(block->log, logInfo, "ResidualScalingUpdated", "", block->label); if (block->callbacks->log_options.log_level >= 5) { jmi_log_node_t inner = jmi_log_enter_vector_(block->log, outer, logInfo, "scaling"); realtype* res = scale_ptr; for (i=0;ilog, 1/res[i]); jmi_log_leave(block->log, inner); } jmi_log_leave(block->log, outer); } } else if(bsop->residual_equation_scaling_mode == jmi_residual_scaling_manual) { for(i = 0; i < N; i++) { scale_ptr[i] = 1/scale_ptr[i]; } /* Scaling is not updated */ if (block->callbacks->log_options.log_level >= 4) { jmi_log_node_t outer = jmi_log_enter_fmt(block->log, logInfo, "ResidualScalingUpdated", "", block->label); if (block->callbacks->log_options.log_level >= 5) { jmi_log_node_t inner = jmi_log_enter_vector_(block->log, outer, logInfo, "scaling"); realtype* res = scale_ptr; for (i=0;ilog, 1/res[i]); jmi_log_leave(block->log, inner); } jmi_log_leave(block->log, outer); } } return; } void jmi_setup_f_residual_scaling(jmi_block_solver_t *block) { jmi_block_solver_options_t* bsop = block->options; int i, N = block->n; realtype* dummy = 0; /* Read manual scaling from annotations and propagated nominal values for * residuals and put them in residual_nominal & scale_ptr*/ if (bsop->residual_equation_scaling_mode == jmi_residual_scaling_manual || bsop->residual_equation_scaling_mode == jmi_residual_scaling_hybrid) { block->F(block->problem_data,dummy, block->residual_nominal, JMI_BLOCK_EQUATION_NOMINAL); for (i = 0; i < N; i++) { if(block->residual_nominal[i] != 0.0) { if(block->residual_nominal[i] < 1/bsop->max_residual_scaling_factor) { block->residual_nominal[i] = 1/bsop->max_residual_scaling_factor; jmi_log_node(block->log, logWarning, "MaxScalingUsed", "Using maximum scaling factor in , " " since specified residual nominal is too small.", block->label, i); } else if(block->residual_nominal[i] > 1/bsop->min_residual_scaling_factor) { block->residual_nominal[i] = 1/bsop->min_residual_scaling_factor; jmi_log_node(block->log, logWarning, "MinScalingUsed", "Using minimal scaling factor in , " " since specified residual nominal is too big.", block->label, i); } } else if(bsop->residual_equation_scaling_mode == jmi_residual_scaling_manual) { block->residual_nominal[i] = 1.0; } } } block->F(block->problem_data,dummy, block->residual_heuristic_nominal, JMI_BLOCK_EQUATION_NOMINAL_AUTO); block->F(block->problem_data,dummy, block->residual_heuristic_nominal, JMI_BLOCK_EQUATION_NOMINAL); for (i = 0; i < N; i++) { if (block->residual_heuristic_nominal[i] == 0.0) { /* TODO: This information could be used to see that the system is structually singular */ block->residual_heuristic_nominal[i] = 1.0; } } } double* jmi_solver_get_f_scales(jmi_block_solver_t* block) { return N_VGetArrayPointer(block->f_scale); } double calculate_scaled_residual_norm(double* residual, double *f_scale, int n) { double scaledRes, norm = 0; int i; for(i=0; i= RAbs(scaledRes)?norm:RAbs(scaledRes); } return norm; }
• ## trunk/RuntimeLibrary/src/jmi/jmi_block_solver.h

 r8507 #include "jmi_log.h" #include #include #include #include #include #include #include /** \brief Evaluation modes for the residual function.*/ void jmi_block_solver_init_default_options(jmi_block_solver_options_t* op); /** \brief Update function scaling based on Jacobian information */ void jmi_update_f_scale(jmi_block_solver_t *block); /**< \brief Retrieve residual scales used in solver */ double* jmi_solver_get_f_scales(jmi_block_solver_t* block); /**< \brief Setup residual scaling */ void jmi_setup_f_residual_scaling(jmi_block_solver_t *block); /** \brief Calculate scaled residual max norm */ double calculate_scaled_residual_norm(double* residual, double *f_scale, int n); /** \brief Check and log illegal iv inputs */ int jmi_check_and_log_illegal_iv_input(jmi_block_solver_t* block, double* ivs, int N);
• ## trunk/RuntimeLibrary/src/jmi/jmi_block_solver_impl.h

 r8313 jmi_string_t label; N_Vector f_scale;          /**< \brief Work vector for scaling of f */ realtype scale_update_time; /**< \brief The last time when f scale was updated */ int n;                         /**< \brief The number of iteration variables */ jmi_real_t* x;                 /**< \brief Work vector for the real iteration variables */ jmi_real_t* last_accepted_x;   /**< \brief Work vector for the real iteration variables holding the last accepted vales by the integrator */ DlsMat J;                       /**< \brief The Jacobian matrix  */ DlsMat J_scale;                 /**< \brief Jacobian matrix scaled with xnorm for used for fnorm calculation */ int using_max_min_scaling_flag; /**< \brief A flag indicating if either the maximum scaling is used of the minimum */ jmi_real_t* dx;                /**< \brief Work vector for the seed vector */ jmi_real_t* max;               /**< \brief Max values for iteration variables */ jmi_real_t* nominal;           /**< \brief Nominal values for iteration variables */ jmi_real_t* residual_nominal;   /**< \brief Nominals values for residual variables */ jmi_real_t* residual_heuristic_nominal;   /**< \brief Heuristic nominals values for residual variables */ jmi_real_t* initial;           /**< \brief Initial values for iteration variables */ jmi_real_t* start_set;         /**< \brief If the start value is specified for the iteration variables */ int jacobian_variability;      /**< \brief Variability of Jacobian coefficients: JMI_CONSTANT_VARIABILITY JMI_PARAMETER_VARIABILITY, JMI_DISCRETE_VARIABILITY, JMI_CONTINUOUS_VARIABILITY */