NAME wnnlp -- constrained non-linear optimization package SYNOPSIS #include "wnnlp.h" #define WN_EQ_COMPARISON 1 #define WN_GT_COMPARISON 2 #define WN_LT_COMPARISON 3 typedef struct wn_linear_constraint_type_struct { ... int size; int *vars; int comparison_type; double *weights; double rhs; } *wn_linear_constraint_type; typedef struct wn_nonlinear_constraint_type_struct { ... int size; int *vars; int comparison_type; double offset; ptr client_data; double (*pfunction)(int size,double *values,ptr client_data); void (*pgradient)(double *grad,int size,double *values,ptr client_data); } *wn_nonlinear_constraint_type; void wn_nlp_conj_method ( int &code, double &val_min, double solution_vect[], double delta_vect[], wn_nonlinear_constraint_type objective, wn_sll constraint_list, int num_vars, int conj_iterations, int offset_iterations, double offset_adjust_rate ) void wn_make_linear_constraint ( wn_linear_constraint_type &constraint, int size,double rhs,int comparison_type ) void wn_make_nonlinear_constraint ( wn_nonlinear_constraint_type &constraint, int size,int comparison_type ) extern int wn_nlp_verbose; DESCRIPTION This package is useful for minimizing general functions of many variables subject to general constraints. The function to minimize and the constraints can be linear or non-linear. If they are non-linear, they must be continuous and differentiable. The constraints can be equality or inequality constraints. The optimization algorithm works roughly as follows: the constraints and the objective function are combined to create a master unconstrained objective function. The constraints are handled by adding a quadratic penalty function to the master objective function for each constraint that is violated. The master objective function is optimized using the conjugate-gradient hill-climbing algorithm (see wnconj). The resulting solution probably violates some of the constraints. Offsets are added to the violated constraints in order to push the solution toward compliance with the violated constraints. The conjugate-gradient procedure is repeated. This entire cycle is repeated until convergence. If all goes well, the algorithm converges to a local minimum; this local minimum is not necessarily a global minimum or even a good local minimum. If either the objective function or the constraint functions are not differentiable, the algorithm often gets stuck in a very sub-optimal solution, and no warning is given that this has happened. Proper scaling of the objective function and of the constraint functions is important to ensure fast and accurate convergence of the algorithm. The functions and variables should be scaled so that the partial derivatives of all functions with respect to all variables in the region of the solution is near +1 and -1. Partials which are always very large or always very small should be scaled so as to be in this region. "wn_nlp_conj_method" solves a constrained non-linear minimization problem. "objective" specifies the linear or non-linear objective function to minimize. "constraint_list" is a singly-linked list of linear and non-linear constraints. "objective" and "constraint_list" together comprise the problem to be solved. They must be composed of objects of types "wn_linear_constraint_type" and "wn_nonlinear_constraint_type". "num_vars" is the number of variables in the problem. The "objective" can be of type wn_linear_constraint_type - use a cast to avoid lint warnings. "solution_vect" is a vector which the caller should initialize to contain a vector which is as close a possible to the solution; after termination "solution_vect" contains the solution. "val_min" is set to the value of the objective function of "solution_vect". "code" contains the return status of the run. See section "Diagnostics" for details on the potential values of "code". "conj_iterations" specifies the number of iterations in the conjugate-gradient algorithm. Making it too small results in erratic behavior and non-convergence; making it too big wastes CPU time. Often setting "conj_iterations" to the same value as "len" works well. "offset_iterations" specifies the number of iterations to adjust constraint offsets; often a number like 10 or 20 works well. "offset_adjust_rate" specifies how aggressively constraint offsets are adjusted. 1 is the most agressive, 0 means no adjustment at all. "offset_adjust_rate" should always be in the range (0,1]. For linear problems, or problems which are only "mildly" nonlinear, 1 usually works well. Using a value too large results in erratic non-convergence; using a value too small wastes CPU time. "delta_vect" is a vector of delta values for computing numerical gradients if the caller does not provide symbolic gradients for all of his non-linear constraints. If all non-linear constraints have symbolic gradients, set "delta_vect" to NULL. The values assigned to "delta_vect" should be large enough so that the functions of the numerical difference are distinguishable within the numerical precision of your machine's double type, but otherwise as small as possible in order to minimize error in the derivative caused by higher-order Taylor expansion terms. "wn_make_linear_constraint" makes a linear constraint or objective of the form: weights[0]*solution_vector[vars[0]] + weights[1]*solution_vector[vars[1]] + .... weights[size-1]*solution_vector[vars[size-1]] comparison to "rhs" "size" specifies the number of variables involved in the constraint; "rhs" is the right hand side of the constraint. "rhs" is ignored if the constraint is used as the objective function. "comparison_type" must be one of WN_EQ_COMPARISON, WN_GT_COMPARISON, or WN_LT_COMPARISON. In the case that an objective function is created, the value of "comparison_type" must be WN_EQ_COMPARISON. After creating the linear constraint, the caller must set the vectors "weights" and "vars" of the linear constraint. "vars" is an array of integers of size "size". Each entry specifies the index of a variable which is included in the constraint. "weights" is an array of doubles of size "size". Each entry specifies a coefficient for the corresponding variable of "vars". "wn_make_nonlinear_constraint" makes a non-linear constraint or objective. "size" specifies the number of variables involved in the constraint; "comparison_type" must be one of WN_EQ_COMPARISON, WN_GT_COMPARISON, or WN_LT_COMPARISON. In the case that an objective function is created the value of "comparison_type" must be WN_EQ_COMPARISON. After creating the non-linear constraint, the caller must set "vars", "pfunction", and possibly "pgradient" and "client_data". "vars" is an array of integers of size "size". Each entry specifies the index of a variable which is included in the constraint. "pfunction" is a pointer to a non-linear function of the variables in "vars". The constraint implemented is one of (*pfunction)(...) >= 0 (*pfunction)(...) <= 0 (*pfunction)(...) == 0 depending on the value of "comparison_type". "pgradient" places the gradient of "pfunction" in "grad". If you do not want to bother figuring out the gradient of "pfunction", do not set "pgradient". The gradient will be computed numerically by calling "pfunction" with domain values differing by amounts from "delta_vect". This can be slow. "client_data" is a pointer to user-defined data to be associated with the constraint. "wn_nlp_verbose" controls the amount of status information output by "wn_nlp_conj_method". 0 produces no output; 1 produces some output; 2 produces more; 3 produces the most. CONVERGENCE ISSUES If you are observing that the solution found by the optimizer is outside of constrained space, then it may be necessary to modify the objective function. The basic idea is to modify the objective function so that it does not go -infinity outside of the contrainted area. This modification allows the quadradic penalty functions (added for each constraint) to take effect and thus makes it easier for the optimizer find a solution within the constrained region. The modified objective function must satisfy the following properties: 1) The modified objective function does not go to -infinity outside the constrained region. 2) The modified objective function is continuous at the boundary of the constrained region. 3) The slope of the modified objective is continuous at the boundary of the constrained region. The single exception to this rule occurs when the objective function has value +infinity at the boundary of the constrained region. For example suppose that F(x) is the objective function with x > 0 as a constraint and you want to create a new Fnew(x) that converges better than F(x). Let s0 be the value of the slope of F(x) (dF/dx) at x=0. Then defining Fnew(x) = (x > 0.0) ? F(x) : ((x > -1.0) ? (s0*x) : (-s0)) may improve your convergence problems. With complex objective functions that use multiplication, division, log(), pow(), an analysis of the behavior of these operators and functions at the constraint boundaries may be necessary to improve convergence. Below are some simple examples of how these operators and functions can cause convergence problems and modifications to them that can help avoid convergence problems. RECIPRICAL OPERATOR EXAMPLE Let the objective function be F(x, ....) = ... + 1/x and the constraint be x > 0 The convergence issue with the "1/x" term in F is that as x approaches 0 from the negative side, F goes to -infinity. Observe that as x approaches 0 from the positive side, F goes to +infinity. Therefore, we should replace 1/x by this function double my_reciprical(double x) { if(x > 0.0) { /* inside constrained region */ return(1.0/x); } else { /* outside constrained region: match function value at boundary of constrained region */ return(FLT_MAX); } } MULTIPLICATION OPERATOR EXAMPLE Let the objective function be F(x, y, ....) = .... + x*y and the constraints be x > 0 y > 0 The convergence issues are: 1) When x > 0 and y approaches -infinity, F goes to -infinity 2) When y > 0 and x approaches -infinity, F goes to -infinity To prevent the multiply from going to -infinity and preserve function values and slopes at the boundary, we replace the multiply with this function double my_multiply(double x,double y) { double prod; prod = x*y; /* This if statement takes effect outside the constrained region only, and so leaves objective function value and slope at boundary of constrained region unchanged */ if(!(prod >= -1.0)) { prod = -1.0; } return(prod); } DIVISION OPERATOR EXAMPLE Let the objective function be F(x, y, ...) = .... + x/y and the constraints be x > 0 y > 0 There are two convergence issues with the "x/y" term in F: 1) As y approaches 0 from the positive side, then F goes to +infinity (if x > 0) 2) When y > 0 and x approaches -infinity F goes to -infinity Note that x/y = x*(1/y) So we can use the above examples to replace the division with double my_divide(double x,double y) { return(my_multiply(x,my_reciprical(y))) } LOG FUNCTION EXAMPLE Let the objective function be F(x,....) = .... + log(x) With constraint x > 1 The convergence issues are: 1) The log function is not defined for x <= 0 2) When x approaches 0 for the positive side, F goes to -infinity Observe that at the boundary of the constraint x = 1, log(x) = 0 and d log(x)/d x = 1. So we can create a better converging F by replacing log(x) by this function double my_log(double x) { if(x > 1.0) { return(log(x)) } else if(x > 0.0) { /* match function value and slope at boundary of constrained region */ return(1.0*(x-1.0)); } else { /* match function value at x = 0, but stop going to -infinity */ return(-1.0); } } POW FUNCTION EXAMPLE Let the objective function be F(x,y,....) = .... + pow(x,y) With constraints x > 0 y > 0 The convergence issues are: 1. When x<0, F becomes imaginary (or rather, undefined in C language) 2. When y<0 and x-->0, F goes to -infinity or +infinity or is undefined Therefore, we replace pow(...) by this function double my_pow(double x,double y) { if(x > 0.0) { return(pow(x,y)); } else { /* match function value at boundary of constrained region. It is impossible to match slope because slope = +infinity at x=0 */ return(0.0); } } RESOURCES wn_nlp_conj_method runs with time = conj_iterations*offset_iterations* stack memory = 1 dynamic memory = len Of course, this figure tells one nothing about the time required to achieve good convergence. This depends heavily on the problem being solved and the user's skill in properly adjusting "conj_iterations", "offset_iterations", and "offset_adjust_rate", and the problem scaling. Typically, if the problem is reasonably behaved and well scaled, one can set "conj_iterations" to something of order of the number of variables "len", and set "offset_iterations" to some number like 10 or 20. The author has succesfully used this routine to solve non-linear problems with 10,000 variables and 50,000 constraints in a few hours of workstation time. This routine can also solve large linear programming problems, but the numerical accuracy is not as good as that of a more specialized LP algorithm. DIAGNOSTICS code == WN_SUCCESS means successful completion, optimum found code == WN_SUBOPTIMAL means termination before optimum found code == WN_UNBOUNDED means function is unbounded EXAMPLES The file "wnlib/acc/conjdir/examples.c.nlp uses the "wnnlp" non-linear constrained optimization package to solve the minimum area buffering problem. The problem formulation is below. PROBLEM FORMULATION: CHOOSE w[i] (1 <= i <= num_vars) TO MINIMIZE sum_over(i){ w[i] } WHILE SATISFYING THE CONSTRAINTS for_all(i) { w[i] >= 0 } sum_over(0<=i<=NUM_VARIABLE_BUFFERS){ w[i+1]/w[i] } <= MAX_DELAY w[0] == DRIVER_BUFFER_SIZE w[NUM_VARIABLE_BUFFERS+1] == RECEIVER_BUFFER_SIZE BUGS SEE ALSO wnconj, wnsplx, wnsll AUTHOR Will Naylor