{VERSION 4 0 "IBM INTEL NT" "4.0" } {USTYLETAB {CSTYLE "Maple Input" -1 0 "Courier" 0 1 255 0 0 1 0 1 0 0 1 0 0 0 0 1 }{CSTYLE "2D Math" -1 2 "Times" 0 1 0 0 0 0 0 0 2 0 0 0 0 0 0 1 }{CSTYLE "2D Output" 2 20 "" 0 1 0 0 255 1 0 0 0 0 0 0 0 0 0 1 } {CSTYLE "" -1 256 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 257 "Courier" 1 10 0 0 0 0 0 1 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 258 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 259 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 260 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 261 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 262 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 263 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{PSTYLE "Normal" -1 0 1 {CSTYLE "" -1 -1 "Tim es" 1 12 0 0 0 1 2 2 2 2 2 2 1 1 1 1 }1 1 0 0 0 0 1 0 1 0 2 2 0 1 } {PSTYLE "Heading 1" -1 3 1 {CSTYLE "" -1 -1 "Times" 1 18 0 0 0 1 2 1 2 2 2 2 1 1 1 1 }1 1 0 0 8 4 1 0 1 0 2 2 0 1 }{PSTYLE "Maple Output" -1 11 1 {CSTYLE "" -1 -1 "Times" 1 12 0 0 0 1 2 2 2 2 2 2 1 1 1 1 }3 3 0 0 0 0 1 0 1 0 2 2 0 1 }{PSTYLE "Author" -1 19 1 {CSTYLE "" -1 -1 " Times" 1 12 0 0 0 1 2 2 2 2 2 2 1 1 1 1 }3 1 0 0 8 8 1 0 1 0 2 2 0 1 } {PSTYLE "Heading 1" -1 256 1 {CSTYLE "" -1 -1 "Times" 1 18 0 0 0 1 2 1 2 2 2 2 1 1 1 1 }3 1 0 0 8 4 1 0 1 0 2 2 0 1 }} {SECT 0 {EXCHG {PARA 256 "" 0 "" {TEXT -1 4 "The " }{TEXT 256 20 "Nonl inearProgramming" }{TEXT -1 8 " Package" }}{PARA 19 "" 0 "" {TEXT -1 58 "by Jason Schattman, Waterloo Maple, Inc., October 16, 2000" }} {PARA 19 "" 0 "" {TEXT -1 24 "jschattman@maplesoft.com" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 8 "restart;" }}}{SECT 0 {PARA 3 "" 0 " " {TEXT -1 16 "The Package Code" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 30396 "NonlinearProgramming := module()\n\nexport\011Optimize, Pri malDualLogBarrier, UnconstrainedNewton;\n\nlocal\011QuadraticIsConvex, FunctionIsQuadratic, FunctionIsLinear, SetIsLinear, \n\011IsPositiveV ector, zerorhs, EvaluateVector, EvaluateDerivative, EvaluateDerivative Sum, \n\011VectorSubtractFast, MatrixSubtractFast, VectorAddFast, Matr ixAddFast, ScaledVectorAdd, \n\011ScaledMatrixAdd, MultiplyMM, Multipl yMV, MultiplyMS, printConstraints, \n\011MultiplyVS, LinearSolveFast, \+ TransposeFast, MatrixColumn, VectorNormFast, DiagMatrix, \n\011SubVect orFast, DotProd, IsPosDefinite, DimensionFast, ZeroVectorFast, ScaledV ectorAddInPlace, \011ScaledMatrixAddInPlace, AddBoundConstraints, GetC onstraintFunctions, ConvertConstraints, PickAlgorithm, OptimizeSimplex , OptimizeExplicitFONOC,\n\011OptimizeInteriorPointMethod, ShouldStop, InitialBarrierParameter, \n\011ComputeBarrierParameter, ComputeDirect ionXW, ComputeDirectionY, \n\011ComputeDirectionW, FeasibleStepSize, M eritStepSize, ConstrainedArmijoStepSize, \n\011ComplexArmijoStepSize, \+ SetMeritPenaltyParameter, ShouldTerminate;\n\noption package;\n\n\nQua draticIsConvex := proc( hess )\n option inline;\n LinearAlgebra:-LA_ Main:-IsDefinite( hess, 'query=positive_def' );\nend proc:\n\nIsPosDef inite := proc( M )\n option inline;\n LinearAlgebra:-LA_Main:-IsDefi nite( M, 'query=positive_def' );\nend proc:\n\nFunctionIsQuadratic := \+ proc( f, n )\n option inline;\n type( f, quadratic( [seq( x[i], i=1. .n )] ) );\nend proc:\n\nFunctionIsLinear := proc( f, n )\n option in line;\n type( f, linear( [seq( x[i], i=1..n )] ) );\nend proc:\n\nSet IsLinear := proc( functionSet, n )\n local i;\n\n for i from 1 to no ps( functionSet ) do\n if not FunctionIsLinear( functionSet[i], n ) then\n return false;\n end if;\n end do;\n\n true;\nend pro c:\n\nIsPositiveVector := proc( v, dim )\n local i;\n\n for i from 1 to dim do\n if v[i] < 0 then\n return false;\n end if;\n \+ end do;\n\n true;\nend proc: \n \nzerorhs := proc( inequal ) # \+ converts the inequality to form h( x )>=0 and returns h( x )\n option inline;\n -`-`( op( inequal ) );\nend proc:\n\nEvaluateVector := pro c( vector_expr, arguments, dim )\n option inline;\n subs( \{seq( x[i ]=arguments[i], i=1..dim )\}, vector_expr );\nend proc:\n\nEvaluateDer ivative := proc( derivative_expr, arguments, dim )\n local expression ;\n\n if type( derivative_expr, 'Vector' ) then\n expression := Tr ansposeFast( derivative_expr );\n else\n expression := derivative_ expr;\n end if;\n\n subs( \{seq( x[i]=arguments[i], i=1..dim )\}, ex pression );\nend proc:\n\nEvaluateDerivativeSum := proc( derivative_ex pr, arguments1, arguments2, dim1, dim2 )\n local expression;\n\n if \+ type( derivative_expr, 'Vector' ) then\n expression := TransposeFas t( derivative_expr );\n else\n expression := derivative_expr;\n e nd if;\n\n subs( \{seq( x[i]=arguments1[i], i=1..dim1 ), seq( k[i]=ar guments2[i], i=1..dim2 )\}, \n expression );\nend proc:\n\nVect orSubtractFast := proc( v1, v2 )\n option inline;\n LinearAlgebra:-L A_Main:-VectorAdd( v1, v2, 1, -1, 'inplace=false', 'outputoptions'=[] \+ );\nend proc:\n\nMatrixSubtractFast := proc( M1, M2 )\n option inline ;\n LinearAlgebra:-LA_Main:-MatrixAdd( M1, M2, 1, -1, 'inplace=false' , 'outputoptions'=[] );\nend proc:\n\nVectorAddFast := proc( v1, v2 ) \n option inline;\n LinearAlgebra:-LA_Main:-VectorAdd( v1, v2, 1, 1, 'inplace=false', 'outputoptions'=[] );\nend proc:\n\nMatrixAddFast := proc( M1, M2 )\n option inline;\n LinearAlgebra:-LA_Main:-MatrixAdd ( M1, M2, 1, 1, 'inplace=false', 'outputoptions'=[] );\nend proc:\n\nS caledVectorAdd := proc( v1, v2, scale1, scale2 )\n option inline;\n \+ LinearAlgebra:-LA_Main:-VectorAdd( v1, v2, scale1, scale2, 'inplace=fa lse', \n 'outputoptions'=[] );\nen d proc:\n\nScaledMatrixAdd := proc( M1, M2, scale1, scale2 )\n option inline;\n LinearAlgebra:-LA_Main:-MatrixAdd( M1, M2, scale1, scale2, 'inplace=false', \n 'outputoption s'=[] );\nend proc:\n\nMultiplyMM := proc( M1, M2 )\n option inline; \n LinearAlgebra:-LA_Main:-MatrixMatrixMultiply( M1, M2, 'inplace=fal se', \n 'outputoptions' =[] );\nend proc:\n\nMultiplyMV := proc( M1, M2 )\n option inline;\n \+ LinearAlgebra:-LA_Main:-MatrixVectorMultiply( M1, M2, 'outputoptions' =[] );\nend proc:\n\nMultiplyMS := proc( M1, s )\n option inline;\n \+ LinearAlgebra:-LA_Main:-MatrixScalarMultiply( M1, s, 'inplace=false', \+ 'outputoptions'=[] );\nend proc:\n\nMultiplyVS := proc( v1, s )\n opt ion inline;\n LinearAlgebra:-LA_Main:-VectorScalarMultiply( v1, s, 'i nplace=false', 'outputoptions'=[] );\nend proc:\n\nLinearSolveFast := \+ proc( M, v )\n option inline;\n LinearAlgebra:-LA_Main:-LinearSolve( M, v, 'method=Cholesky', 'free=s', 'conjugate=false', \n \+ 'inplace=false', 'outputoptions'=[] );\nend pr oc:\n\nTransposeFast := proc( M )\n option inline;\n LinearAlgebra:- LA_Main:-Transpose( M, 'inplace=false', 'outputoptions'=[] );\nend pro c:\n\nMatrixColumn := proc( M, c )\n option inline;\n LinearAlgebra: -LA_Main:-Column( M, c, 'outputoptions'=[] );\nend proc:\n\nVectorNorm Fast := proc( v )\n option inline;\n LinearAlgebra:-LA_Main:-VectorN orm( v, 2 );\nend proc:\n\nDiagMatrix := proc( veclist, n)\n option i nline;\n LinearAlgebra:-LA_Main:-DiagonalMatrix( veclist, n, n,'outpu toptions'=[datatype=float[8]] );\nend proc:\n\nSubVectorFast := proc( \+ v, index )\n option inline;\n LinearAlgebra:-LA_Main:-SubVector( v, \+ index, 'outputoptions'=[] );\nend proc:\n\nDotProd := proc( v1, v2 )\n option inline;\n LinearAlgebra:-LA_Main:-DotProduct( v1, v2, 'conju gate=false' );\nend proc:\n\nDimensionFast := proc( M )\n option inli ne;\n LinearAlgebra:-LA_Main:-Dimension( M );\nend proc:\n\n \nZeroV ectorFast := proc( dim )\n option inline;\n LinearAlgebra:-LA_Main:- ZeroVector[column]( dim, 'compact=false', 'outputoptions'=[] );\nend p roc:\n \nScaledVectorAddInPlace := proc( v1, v2, scale1, scale2 )\n \+ option inline;\n LinearAlgebra:-LA_Main:-VectorAdd( v1, v2, scale1, s cale2, 'inplace=true', \n 'outputo ptions'=[] );\nend proc:\n\nScaledMatrixAddInPlace := proc( M1, M2, sc ale1, scale2 )\n option inline;\n LinearAlgebra:-LA_Main:-MatrixAdd( M1, M2, scale1, scale2, 'inplace=true', \n \+ 'outputoptions'=[] );\nend proc:\n\nAddBoundConstraints := p roc( constraintSet::set, variableBounding, numVars::integer )\n \n i f variableBounding = 'free' then\n constraintSet;\n elif variableB ounding = 'nonnegative' then\n constraintSet union \{seq(x[i]>=0, i = 1..numVars)\};\n else\n error \"Variable bounds must have value 'free' or 'nonnegative' but received %2\", \n variableBoundi ng;\n end if;\n\nend proc:\n\n \n \nConvertConstraints := proc( cons traintSet::set, variableBounding, numVars::integer )\n local constrai ntsWithBounds, constraintFunctions;\n \n constraintsWithBounds := Ad dBoundConstraints( constraintSet, variableBounding, numVars );\n cons traintFunctions := GetConstraintFunctions( constraintsWithBounds );\n \n Vector( convert( constraintFunctions, list ) );\nend proc:\n\nGetC onstraintFunctions := proc( constraintSet::set )\n map( zerorhs, cons traintSet );\nend proc:\n\nprintConstraints := proc ( constraints, m ) \n [ seq(constraints[i]>=0, i=1..m) ];\nend proc:\n\nOptimize := proc (function, objective, num_variables, \n constraint_set , variable_bounds, x_start)\n local obj_function, constraint_function s, num_constraints, all_constraints;\n\n if objective='max' then\n \+ obj_function := -1 * function;\n userinfo( 2, 'Optimize', 'Maximiz e' );\n elif objective='min' then\n obj_function := function;\n \+ userinfo( 2, 'Optimize', 'Minimize' );\n else\n error \"expecting the %-1 argument to be either `%2` or `%3`\", \n 2, 'min', 'max' ;\n end if;\n\n userinfo( 1, 'Optimize', function );\n\n PickAlgori thm( obj_function, num_variables, constraint_set, variable_bounds, x_s tart );\nend proc:\n\n\nPickAlgorithm := proc( obj_function, numVars, \+ constraintSet, variableBounding, x_start )\n local n, all_constraints , constraintFunctions, constraintStructure, hess;\n\n n := numVars;\n \n if nops( constraintSet ) = 0 and variableBounding = 'free' then\n \+ constraintStructure := unconstrained;\n else\n constraintStruct ure := constrained;\n end if;\n\n if SetIsLinear( GetConstraintFunct ions( constraintSet ), n ) then\n if FunctionIsLinear( obj_function , n ) then\n OptimizeSimplex( obj_function, constraintSet, variab leBounding, n );\n elif FunctionIsQuadratic( obj_function, n ) then \n hess := Matrix( hessian( obj_function, [seq( x[i], i=1..n )] ) );\n\n if QuadraticIsConvex( hess ) then\n if constraintS tructure = unconstrained then\n OptimizeExplicitFONOC( obj_fu nction, hess, n );\n else\n OptimizeInteriorPointMetho d( obj_function, n, constraintSet, variableBounding, \n \+ x_start, convex, constrained );\n end i f;\n else\n OptimizeInteriorPointMethod( obj_function, n, \+ constraintSet, variableBounding,\n \+ x_start, nonconvex, constraintStructure );\n end if;\n else\n OptimizeInteriorPointMethod( obj_function, n, constraintSet, var iableBounding,\n x_start, nonconvex, constraintStructure );\n end if;\n else\n OptimizeInteriorPoin tMethod( obj_function, n, constraintSet, variableBounding,\n \+ x_start, nonconvex, constraintStructure );\n e nd if;\nend proc:\n\nOptimizeSimplex := proc( objFunction, constraintS et, variableBounding, numVars )\n local constraintsWithBounds;\n\n u serinfo( 2, 'Optimize', `Linear problem: solving with simplex method` \+ );\n\n constraintsWithBounds := AddBoundConstraints( constraintSet, v ariableBounding, numVars ); \n \n simplex['minimize']( objFunction, \+ constraintsWithBounds );\n \nend proc:\n\nOptimizeExplicitFONOC := pr oc( obj_function, hess, n )\n local gradient, gradient0, soln;\n\n u serinfo( 2, 'Optimize', \n `Convex quadratic problem: solvi ng explicitly with first-order conditions` );\n\n gradient := Vector( grad( obj_function, [seq( x[i], i=1..n )] ) );\n gradient0 := Evalua teDerivative( gradient, ZeroVectorFast( n ), n );\n\n userinfo( 2, 'O ptimizeExplicitFONOC', `Global optimum found at` );\n\n soln := Linea rSolveFast( hess, MultiplyVS( gradient0, -1 ) );\n\n \{ seq( x[i] = s oln[i], i = 1..n) \};\nend proc:\n\nOptimizeInteriorPointMethod := pro c( obj_function, numVars, constraintSet, variableBounding, \n \+ x_start, convexityType, constraintStructur e )\n if constraintStructure = constrained then\n userinfo( 2, 'Op timize', \n `Constrained`, convexityType, `problem: solvi ng with`, \n convexityType, `primal-dual log-barrier alg orithm`);\n PrimalDualLogBarrier( obj_function, numVars, constraint Set, variableBounding, \n x_start, convexityT ype);\n\n else\n userinfo( 2, 'Optimize', \n `Unconst rained`, convexityType, `problem: solving with`, \n conv exityType, `Newton's Method`);\n UnconstrainedNewton( obj_function, numVars, x_start, convexityType, float[8] );\n end if;\n\nend proc: \n\nPrimalDualLogBarrier := proc(obj_function, numVars, \n \+ constraintSet, variableBounding, x_start, convexityTy pe)\n local x_curr, w_slack, y_mult, n, m, fx, constraints, \n \+ jacobian_constr_T, jacobian_constr, constraint_start0, constraint_sta rt1, \n hess_constr, hess_obj_fun, hess_lag_fun, grad_obj_fun, \+ grad_lag_fun, \n w_inv, W_inv, Y, W_inv_Y, constraints_curr, gr ad_obj_curr, hess_lag_curr, \n jacobian_curr, grad_lag_curr, rh o, gamma, sigma, y_inv, W_Y_inv, dir_x, \n dir_y, dir_w, b, ste psize_max, stepsize, merit_function, merit, \n grad_merit_fun, \+ B, mu, mu_max, grad_merit_curr, merit_grad_dot_dir, dir_xw, \n \+ infeasibility_tolerance, slack_tolerance, max_barrier_iterations, \n \+ norm_grad_lag_start, norm_w_slack_start, iter_count;\n\n ASSERT ( type( numVars, 'integer' ) );\n ASSERT( nargs = 6 );\n\n constrain ts := ConvertConstraints( constraintSet, variableBounding, numVars ); \n\n n := numVars;\n m := DimensionFast( constraints );\n\n constra int_start0 := subs( \{seq( x[i]=x_start[i], i=1..n )\}, constraints ); \n constraint_start1 := map( x->max( 1, x ), constraint_start0 );\n\n w_slack := Vector( map( evalf, constraint_start1 ), datatype=float[8 ] );\n y_mult := Vector( m, .95, datatype=float[8] );\n x_curr := Ve ctor( map( evalf, x_start ), datatype=float[8] );\n\n userinfo( 1, 'P rimalDualLogBarrier', `Minimize`, print(obj_function) );\n userinfo( \+ 1, 'PrimalDualLogBarrier', `subject to`, printConstraints(constraints, m) );\n userinfo( 2, 'PrimalDualLogBarrier', `Starting point` );\n \+ userinfo( 2, 'PrimalDualLogBarrier', print(x_curr) );\n userinfo( 3, \+ 'PrimalDualLogBarrier', print(Y0), print(y_mult), print(W0), print(w_s lack) );\n\n grad_obj_fun := Vector( grad( obj_function, [seq( x[i], \+ i=1..n )] ) );\n hess_obj_fun := Matrix( hessian( obj_function, [seq( x[i], i=1..n )] ) );\n\n jacobian_constr := Matrix( jacobian( constr aints, [seq( x[i], i=1..n )] ) );\n jacobian_constr_T := TransposeFas t( jacobian_constr );\n hess_constr := map( Matrix, map( hessian, con straints, [seq( x[i], i=1..n )] ) );\n\n grad_lag_fun := VectorSubtra ctFast( \n grad_obj_fun, \n Vector ( evalm( add( MultiplyVS( MatrixColumn( jacobian_constr_T, i ), \n \+ k[i] ), \+ \+ i=1..m ) ) ) );\n hess_lag_fun := MatrixS ubtractFast( \n hess_obj_fun, \n M atrix( evalm( add( MultiplyMS( hess_constr[i], k[i] ), i=1..m ) ) ) ); \n\n merit := obj_function - mu*add( log( k[i] ), i=1..m ) + .5*B*add ( (constraints[i]-k[i])^2,\n \+ i=1..m );\n merit_function := unapply( meri t, seq( x[i], i=1..n ), seq( k[i], i=1..m ), mu, B );\n grad_merit_fu n := Vector( grad( merit, [seq( x[i], i=1..n ), seq( k[i], i=1..m )] ) );\n \n userinfo( 3, 'PrimalDualLogBarrier', GRAD_F, print(Transpos eFast( grad_obj_fun ) ));\n userinfo( 3, 'PrimalDualLogBarrier', HESS IAN_F, print(hess_obj_fun ));\n userinfo( 3, 'PrimalDualLogBarrier', \+ JACOBIAN_OF_CONSTRAINTS, print(jacobian_constr ));\n userinfo( 3, 'Pr imalDualLogBarrier', GRAD_LAGRANGIAN, print(TransposeFast( grad_lag_fu n )) );\n userinfo( 3, 'PrimalDualLogBarrier', HESSIAN_LAGRANGIAN, pr int(hess_lag_fun) );\n userinfo( 3, 'PrimalDualLogBarrier', print(mer it) );\n userinfo( 3, 'PrimalDualLogBarrier', print(grad_merit_fun) ) ;\n\n B := .1;\n max_barrier_iterations := min( 200, 50 * n );\n in feasibility_tolerance := .00001;\n slack_tolerance := .00001 * m;\n \+ \n constraints_curr := evalf( EvaluateVector( constraints, x_curr, n ) );\n grad_lag_curr := evalf( EvaluateDerivativeSum( grad_lag_fun, \+ x_curr, y_mult, n, m ) );\n rho := VectorSubtractFast( w_slack, const raints_curr );\n\n norm_grad_lag_start := VectorNormFast( grad_lag_cu rr );\n norm_w_slack_start := VectorNormFast( w_slack );\n\n mu := I nitialBarrierParameter( w_slack, constraints_curr, n, m );\n mu_max : = 100000 / m;\n\n iter_count := 1;\n\n to max_barrier_iterations \n \+ while not ShouldStop( grad_lag_curr, norm_grad_lag_start, \n \+ rho, norm_w_slack_start, w_slack, y_mult, infeasibilit y_tolerance, \n slack_tolerance ) do\n ja cobian_curr := evalf( EvaluateDerivative( jacobian_constr, x_curr, n ) );\n grad_obj_curr := evalf( EvaluateDerivative( grad_obj_fun, \+ x_curr, n ) );\n hess_lag_curr := evalf( EvaluateDerivativeSum( \+ hess_lag_fun,x_curr, y_mult, n, m ) );\n grad_merit_curr := eval f( EvaluateDerivativeSum(grad_merit_fun,x_curr,w_slack,n,m ));\n\n \+ w_inv := map( x->1/x, w_slack );\n y_inv := map( x->1/x, y_mu lt );\n W_inv_Y := DiagMatrix( [seq( w_inv[i] * y_mult[i], i=1.. m )], m );\n W_Y_inv := DiagMatrix( [seq( w_slack[i] * y_inv[i], i=1..m )], m );\n \n gamma := ScaledVectorAdd( w_inv, y_m ult, mu, -1 );\n b := VectorAddFast( MultiplyMV( W_inv_Y, rho ), gamma );\n\n dir_xw := ComputeDirectionXW( grad_lag_curr, hess_ lag_curr, jacobian_curr, \n W_inv_ Y, b,grad_merit_curr, rho, constraints_curr, \n \+ grad_obj_curr, n, convexityType );\n dir_x := SubV ectorFast( dir_xw, [1..n] );\n dir_w := SubVectorFast( dir_xw, [ n+1..n+m] );\n dir_y := ComputeDirectionY( dir_w, W_inv_Y, b, rh o );\n\n merit_grad_dot_dir := evalf( DotProd( grad_merit_curr, \+ dir_xw ) );\n\n if merit_grad_dot_dir >= 0 then\n B := S etMeritPenaltyParameter( dir_x, grad_obj_curr, rho, mu, jacobian_curr, \n w_inv, m );\n end if; \n \+ \n stepsize := MeritStepSize( merit_function, n, m, x_curr, w_slack, y_mult, \n di r_x, dir_y, dir_w, dir_xw, merit_grad_dot_dir, mu, B );\n\n if s tepsize < .00000001 and iter_count >= 30 then\n userinfo( 2, ' PrimalDualLogBarrier', \n \"Local optimum not found. \+ Best solution found (may not be feasible):\" );\n return x_cu rr;\n elif stepsize < .000001 then\n B := 2 * B;\n \+ end if;\n\n ScaledVectorAddInPlace( x_curr, dir_x, 1, stepsize \+ );\n ScaledVectorAddInPlace( y_mult, dir_y, 1, stepsize ) ;\n \+ ScaledVectorAddInPlace( w_slack, dir_w, 1, stepsize ) ;\n\n \+ userinfo( 3, 'PrimalDualLogBarrier', MU, mu, B_IS, B );\n userin fo( 2, 'PrimalDualLogBarrier', print(TransposeFast( x_curr )) );\n \+ userinfo( 3, 'PrimalDualLogBarrier', print(dir_x), print(x_curr) ); \n userinfo( 3, 'PrimalDualLogBarrier', print(Y_DIR), print(dir_ y), print(y_mult) );\n userinfo( 3, 'PrimalDualLogBarrier', prin t(W_DIR), print(dir_w), print(w_slack) );\n \n constraints _curr := evalf( EvaluateVector( constraints, x_curr, n ) );\n gr ad_lag_curr := evalf( EvaluateDerivativeSum( grad_lag_fun, x_curr, y_m ult, n, m ));\n rho := VectorSubtractFast( w_slack, constraints_ curr );\n userinfo( 4, 'PrimalDualLogBarrier', \n \+ H( X ), constraints_curr, W, w_slack, INFEAS, rho );\n\n mu := min( ComputeBarrierParameter( w_slack, y_mult, m, constraints_curr ), mu_max );\n iter_count := iter_count + 1;\n end do;\n\n if co nvexityType = convex then\n userinfo( 2, 'PrimalDualLogBarrier', `G lobal optimum found at` );\n else\n userinfo( 2, 'PrimalDualLogBar rier', `Local optimum found at` );\n end if;\n\n \{ seq( x[i] = x_cu rr[i], i = 1..n) \};\nend proc: \n\nShouldStop := proc(grad_lag_curr, norm_grad_lag_start, rho, norm_w_slack_start, \n w_ slack, y_mult, infeas_tolerance, slack_tolerance )\n local relative_p rimal_infeas, relative_dual_infeas, relative_comp_slack;\n\n relative _dual_infeas := VectorNormFast( grad_lag_curr ) / norm_grad_lag_start; \n \n relative_primal_infeas := VectorNormFast( rho ) / norm_w_slack _start;\n\n relative_comp_slack := DotProd( w_slack, y_mult );\n\n u serinfo( 3, 'PrimalDualLogBarrier', \n REL_DUAL_INFEAS, rel ative_dual_infeas, REL_PRIMAL_INFEAS, \n relative_primal_in feas, REL_COMP_SLACK, relative_comp_slack );\n\n evalb( relative_dual _infeas <= infeas_tolerance and \n relative_primal_infeas <= i nfeas_tolerance and \n relative_comp_slack <= slack_tolerance \+ );\nend proc:\n\nInitialBarrierParameter := proc( w_slack, constraints _curr, n, m )\n local big_M, mu_init;\n\n big_M := 200 / m;\n\n mu_ init := big_M * max( VectorNormFast( w_slack ), VectorNormFast( constr aints_curr ) );\n userinfo( 2, 'Optimize', `Initial barrier parameter :`, mu_init );\n\n mu_init;\nend proc:\n\nComputeBarrierParameter := \+ proc(w_slack, y_mult, m, constraints_curr)\n local dotprod_wy, mu_def ault, wy, min_wy, xi, i, scale;\n\n wy := seq( w_slack[i] * y_mult[i] , i=1..m );\n dotprod_wy := `+`( wy );\n\n scale := .1;\n\n mu_defa ult := scale * dotprod_wy / m;\n\n min_wy := min( wy );\n\n xi := s cale * min_wy / mu_default;\n\n if xi = 1 then\n mu_default;\n el if xi = 0 then\n 8 * mu_default;\n else\n min( 2, .95*( 1-xi )/ xi )^3 * mu_default;\n end if;\nend proc:\n\nComputeDirectionXW := pr oc(grad_lag_curr, hess_lag_curr, jacobian_curr, \n \+ W_inv_Y, b, grad_merit_curr, rho, \n co nstraints_curr, grad_obj_curr, n, \n convexity Type )\n local W_inv_YA, A_trans_W_inv_YA, N, v, H, H_diag, lambda, d ir_x, dir_w, i;\n\n W_inv_YA := MultiplyMM( W_inv_Y, jacobian_curr ); \n A_trans_W_inv_YA := MultiplyMM( TransposeFast( jacobian_curr ), W_ inv_YA );\n N := MatrixAddFast( hess_lag_curr, A_trans_W_inv_YA );\n \+ v := VectorSubtractFast( MultiplyMV( TransposeFast( jacobian_curr ), \+ b ), \n grad_lag_curr );\n\n if convexityTy pe = nonconvex then\n H := hess_lag_curr;\n H_diag := DiagMatrix ( [seq( max(abs( H[i, i] ),.01), i=1..n )], n );\n lambda := 1.2;\n i := 1;\n\n to 20 while IsPosDefinite( N ) = false do \n S caledMatrixAddInPlace( N, H_diag, 1, lambda );\n userinfo(3,'Opti mize', \n `Generalized Hessian adjusted by diagonal fact or`, lambda );\n if i < 10 then\n lambda := 2 * la mbda;\n else\n lambda := 10 * lambda;\n e nd if;\n end do;\n end if;\n\n dir_x := LinearSolveFast( N, v ); \n\n dir_w := VectorSubtractFast( MultiplyMV( jacobian_curr, dir_x ), rho );\n\n Vector( map( evalf, [dir_x, dir_w] ), datatype=float[8] ) ;\nend proc:\n \nComputeDirectionY := proc(dir_w, W_inv_Y, b, rho ) \n local dir_y;\n\n dir_y := VectorSubtractFast( b, MultiplyMV( W_in v_Y, VectorAddFast( dir_w, rho ) ) );\n\n Vector( map( evalf, dir_y ) , datatype=float[8] );\nend proc:\n\nComputeDirectionW := proc(dir_y, \+ W_Y_inv, gamma) #not currently called\n MultiplyMV( W_Y_inv, VectorSu btractFast( gamma, dir_y ) );\nend proc:\n\nFeasibleStepSize := proc(d ir_w, dir_y, w_slack, y_mult, \n m)\n local alp ha_max, w_ratio, y_ratio, max_ratio, min_ratio;\n\n max_ratio := max( seq( -1*dir_w[i]/w_slack[i], i=1..m ), \n seq( -1* dir_y[i]/y_mult[i], i=1..m ) );\n if max_ratio <= 1 then\n 0.95;\n else\n 0.95 / max_ratio;\n end if;\nend proc:\n\nMeritStepSize : = proc(merit_fun, n, m, x_curr, \n w_slack, y_mul t, \n dir_x, dir_y, \n dir_w , dir_xw, merit_grad_dot_dir, mu, B)\n local stepsize_max, stepsize_a rmijo, grad_merit_curr, xw_curr, \n direction_xw, grad_merit_xw , x_new, constraints_new;\n\n stepsize_max := FeasibleStepSize( dir_w , dir_y, w_slack, y_mult, m );\n\n xw_curr := Vector( [x_curr, w_slac k, mu, B] );\n direction_xw := Vector( [dir_xw, 0, 0] );\n\n stepsiz e_armijo := ConstrainedArmijoStepSize( merit_fun, n+m+2, xw_curr, dire ction_xw, \n merit_grad _dot_dir, stepsize_max );\n userinfo( 3, 'NoName', \"Max stepsize %g, Armijo stepsize %g\\n\", stepsize_max, stepsize_armijo );\n\n stepsi ze_armijo;\nend proc:\n\nConstrainedArmijoStepSize := proc(F, dim, poi nt_curr, \n direction, grad_dot_direc tion, stepsize_max)\n local beta, sigma, m, point_new, improvement, t hreshold, beta_m;\n\n beta := .5;\n beta_m := stepsize_max;\n sigma := .001;\n m := 0;\n\n point_new := VectorAddFast( point_curr, dire ction );\n\n improvement := F( seq( point_curr[i], i=1..dim ) ) - F( \+ seq( point_new[i], i=1..dim ) );\n threshold := -1 * sigma * grad_dot _direction;\n\n to 50 while improvement < threshold and abs( improvem ent ) > 0 do\n beta_m := beta * beta_m;\n point_new := ScaledVec torAdd( point_curr, direction, 1, beta_m );\n improvement := F( seq ( point_curr[i], i=1..dim ) ) \n - F( seq( point_new [i], i=1..dim ) );\n threshold := evalf( beta * threshold );\n m := m+1;\n end do;\n\n beta_m;\nend proc:\n \nComplexArmijoStepSize := proc(F, dim, point_curr, \n direction , grad_dot_direction, stepsize_max)\n local beta, sigma, m, point_new , improvement, realImprovement, threshold, beta_m;\n\n beta := .5;\n \+ beta_m := stepsize_max;\n sigma := .001;\n m := 0;\n\n point_new : = VectorAddFast( point_curr, direction );\n\n improvement := evalf( F ( seq( point_curr[i], i=1..dim ) ) )\n - evalf( F( seq ( point_new[i], i=1..dim )) );\n realImprovement := Re( improvement ) ;\n\n threshold := -1 * sigma * grad_dot_direction;\n\n to 50 while \+ realImprovement < threshold and abs( realImprovement ) > 0 do\n bet a_m := beta * beta_m;\n point_new := ScaledVectorAdd( point_curr, d irection, 1, beta_m );\n improvement := evalf( F( seq( point_curr[i ], i=1..dim ) ) )\n - evalf( F( seq( point_new[i], i =1..dim )) );\n realImprovement := Re( improvement );\n threshol d := evalf( beta * threshold );\n m := m+1;\n end do;\n\n beta_m; \nend proc:\n\nSetMeritPenaltyParameter := proc(dir_x, grad_obj_curr, \+ rho, mu, \n jacobian_curr, w_inv, m)\n local v, k, jacob_T, squared_norm_rho;\n\n jacob_T := TransposeFast ( jacobian_curr );\n\n v := ScaledVectorAdd( grad_obj_curr, MultiplyM V( jacob_T, w_inv ), 1, -1*mu );\n k := mu * DotProd( w_inv, rho );\n \n squared_norm_rho := add( rho[i]^2, i=1..m );\n\n 2 * ( DotProd( v , dir_x ) + k ) / squared_norm_rho;\nend proc: \n\n\nUnconstrainedNew ton := proc( f, n, x_start, convexityType, datType )\n local x_curr, \+ direction, grad_fun, hess_fun, grad_start, H, H_diag, lambda, \n \+ norm_grad_start, grad_curr, norm_grad_curr, hess_curr, hess_inv, \+ \n neg_grad_curr, i, fx, f_curr, step_size, grad_d ot_direction, max_step_size, \n max_iterations;\n\n if not (da tType = 'float[8]' or datType = 'complex[8]') then \n error \"dat atype must be float[8] or complex[8] but received `%2`\", 2, datType; \n end if;\n\n fx := unapply( f, seq( x[i], i=1..n ) );\n grad_fun \+ := Vector['row']( grad( f, [seq( x[i], i=1..n )] ) );\n hess_fun := M atrix( hessian( f, [seq( x[i], i=1..n )] ) );\n\n x_curr := Vector( m ap( evalf, x_start ), datatype=float[8] );\n\n grad_start := evalf( m ap( x->Re(x), EvaluateDerivative( grad_fun, x_curr, n )) );\n norm_gr ad_start := VectorNormFast( grad_start );\n\n if norm_grad_start = 0 \+ and convexityType = convex then\n return x_start;\n end if;\n\n u serinfo( 2, 'UnconstrainedNewton', \"Starting point\" );\n userinfo( \+ 2, 'UnconstrainedNewton', print(x_curr) );\n userinfo( 3, 'Unconstrai nedNewton', Y0, y_mult, W0, w_slack );\n userinfo( 3, 'UnconstrainedN ewton', GRAD_F, grad_fun );\n userinfo( 3, 'UnconstrainedNewton', HES SIAN_F, hess_fun );\n \n grad_curr := grad_start;\n direction := Mu ltiplyVS( grad_curr, -1 );\n\n max_iterations := 100;\n\n if datType = 'float[8]' then\n to max_iterations while not ShouldTerminate( n orm_grad_start, grad_curr, direction ) do\n neg_grad_curr := Mult iplyVS( grad_curr, -1 );\n hess_curr := evalf( EvaluateDerivative ( hess_fun, x_curr, n ) );\n\n H := hess_curr;\n\n if convex ityType = nonconvex then\n H_diag := DiagMatrix( [seq( max(abs( H[i, i] ),.01), i=1..n )], n );\n lambda := 1.2;\n\n i \+ := 1;\n to 20 while IsPosDefinite( H ) = false do \n S caledMatrixAddInPlace( H, H_diag, 1, lambda );\n userinfo(3,' UnconstrainedNewton',\n `Hessian increased by diagon al factor`,lambda);\n if i < 10 then\n lambda := 2 * lambda;\n else\n lambda := 10 * lambda;\n \+ end if;\n i := i + 1;\n end do;\n end if;\n\n direction := LinearSolveFast( H, neg_grad_curr );\n\n grad_ dot_direction := DotProd( grad_curr, direction );\n \n max_step _size := .95;\n \n step_size := ConstrainedArmijoStepSize( fx, n, x_curr, direction, grad_dot_direction, \n \+ max_step_size );\n\n ScaledVectorAddInPlace( x_curr, direction, 1, step_size );\n grad_curr := EvaluateDerivative( gr ad_fun, x_curr, n );\n\n userinfo(2, 'UnconstrainedNewton', print (TransposeFast( x_curr )) );\n userinfo(3, 'UnconstrainedNewton', X_DIR, print(direction), STEPSIZE, step_size, X, print(x_curr) );\n \+ end do:\n\n else # FOR COMPLEX MATRICES\n to max_iterations whil e not ShouldTerminate( norm_grad_start, grad_curr, direction ) do\n \+ neg_grad_curr := evalf( MultiplyVS( map( x->Re(x), grad_curr), -1 ) );\n\n hess_curr := evalf( map( x->Re(x), EvaluateDerivative( he ss_fun, x_curr, n ) ) );\n\n H := hess_curr;\n\n if convexit yType = nonconvex then\n H_diag := DiagMatrix( [seq( max(abs( H [i, i] ),.01), i=1..n )], n );\n lambda := 1.2;\n\n i := 1;\n to 20 while IsPosDefinite( H ) = false do \n Sca ledMatrixAddInPlace( H, H_diag, 1, lambda );\n userinfo(3,'Un constrainedNewton',\n `Hessian increased by diagonal factor`, lambda);\n if i < 10 then\n lambda := 2 \+ * lambda;\n else\n lambda := 10 * lambda;\n \+ end if;\n i := i + 1;\n end do;\n end if;\n\n \+ direction := LinearSolveFast( H, neg_grad_curr ); \n grad_do t_direction := DotProd( grad_curr, direction );\n\n max_step_size := .95;\n \n step_size := ComplexArmijoStepSize( fx, n, x_curr, \+ direction, \n grad_dot_direct ion, max_step_size );\n\n ScaledVectorAddInPlace( x_curr, directi on, 1, step_size );\n grad_curr := evalf( EvaluateDerivative( gra d_fun, x_curr, n ) );\n\n userinfo(2, 'UnconstrainedNewton', prin t(TransposeFast( x_curr )) );\n userinfo(3, 'UnconstrainedNewton' , X_DIR, print(direction), STEPSIZE,step_size,X, print(x_curr) );\n \+ end do:\n\n end if;\n\n if convexityType = convex then\n userinf o( 2, 'UnconstrainedNewton', `Global optimum found at` );\n else\n \+ userinfo( 2, 'UnconstrainedNewton', `Local optimum found at` );\n en d if;\n\n \{ seq( x[i] = x_curr[i], i = 1..n) \};\nend proc:\n\nShoul dTerminate := proc( norm_grad_start, grad_curr, direction )\n local n orm_grad_curr, norm_direction, rel_grad_norm;\n\n norm_grad_curr := e valf( VectorNormFast( grad_curr ) );\n norm_direction := evalf( Vecto rNormFast( direction ) );\n\n rel_grad_norm := norm_grad_curr / norm_ grad_start;\n\n userinfo( 3, 'UnconstrainedNewton', \n REL _GRAD_NORM, print(rel_grad_norm), \n REL_DIR_X, print(norm_ direction) );\n\n evalb( rel_grad_norm < .001 and norm_direction < .0 01 );\nend proc:\n\nend module:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}{SECT 0 {PARA 3 "" 0 "" {TEXT -1 11 "Saving the " }{TEXT 258 21 "NonlinearProgramming " }{TEXT -1 22 "Package into a Library" } }{EXCHG {PARA 0 "" 0 "" {TEXT -1 12 "To save the " }{TEXT 260 20 "Nonl inearProgramming" }{TEXT -1 65 " package into a library, remove the sh arp character (#) from the " }{TEXT 261 7 "libname" }{TEXT -1 22 " ass ignment below and " }{TEXT 259 0 "" }{TEXT -1 154 "replace the example directory name \"C:/mylib/nlp\" with the name of the directory where \+ you want the package saved. Then remove the # characters from the " } {TEXT 262 5 "march" }{TEXT -1 5 " and " }{TEXT 263 8 "savelib " } {TEXT -1 17 "commands. \n " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 33 "libname:=\"C:/mylib/nlp\",libname; " }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%(libnameG6$Q-C:/mylib/nlp6\"Q=C:\\Program~Files\\Mapl e~6/libF'" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 32 "#march('create ',libname[1],100);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 32 "savel ib('NonlinearProgramming');" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 29 "with( NonlinearProgramming );" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#7 %%)OptimizeG%5PrimalDualLogBarrierG%4UnconstrainedNewtonG" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 57 "\nNow execute the worksheet by choosing from the menu bar " }} {PARA 0 "" 0 "" {TEXT 257 24 "Edit->Execute->Worksheet" }}}}}{MARK "2 \+ 1 0 0" 20531 }{VIEWOPTS 1 1 0 1 1 1803 1 1 1 1 }{PAGENUMBERS 0 1 2 33 1 1 }