## ams_version=1.0 Model Main_Multi_Start { Comment: { "Keywords: Nonlinear Programming, NLP, Multistart Algorithm, GMP, Module." } Section The_Model { DeclarationSection Actual_Model_Declaration { Parameter ARandomNumber { Text: "Changed when regenerate the function"; InitialData: 5; } Variable x { Range: [0, 5]; } Variable y { Definition: sin(ARandomNumber*x) * sin(15*x) * (x - 2.5); } MathematicalProgram MP { Objective: y; Direction: maximize; } } DeclarationSection Multi_Start_Settings_Declaration { Parameter MaxIteration { Text: "Limits the number of iterations by the MultiStart algorithm"; Range: { {1..inf} } InitialData: 25; } Parameter SamplePoints { Text: "Number of sample points to calculate penalized objective"; Range: integer; InitialData: 25; } Parameter SelectedSamplePoints { Text: "Number of sample points selected as start points to solve the NLP"; Range: { {1..SamplePoints} } InitialData: 10; } Parameter ShowProgress { Text: "Show Multistart Progress"; InitialData: 0; } Parameter ShowSolvedPoints { Text: "Show Solved Points"; InitialData: 1; } Parameter ShowAssignedPoints { Text: "Show Assigned Points"; InitialData: 0; } } DeclarationSection MultStart_Results_Declaration { Parameter AllSolvesCounter; Set AllSolves { SubsetOf: Integers; Index: as, as2; } Parameter StartPoint { IndexDomain: (as); } Parameter EndPoint { IndexDomain: (as); Definition: sin(5*StartPoint(as)) * sin(15*StartPoint(as)) * (StartPoint(as) - 2.5); } Parameter EndPoint_x { IndexDomain: (as); } Parameter EndPoint_y { IndexDomain: (as); } } } Section GUI { DeclarationSection Curve_Declarations_Declaration { Set DiscretePoints { Index: xp; Parameter: SelectedPoint; Definition: Data{1..501}; } Parameter DP_x { IndexDomain: (xp); Definition: 5/[card(DiscretePoints)-1] * [ord(xp)-1]; } Parameter DP_y { IndexDomain: (xp); Definition: sin(ARandomNumber*DP_x(xp)) * sin(15*DP_x(xp)) * (DP_x(xp)- 2.5); } Parameter CurveMarkers { IndexDomain: (xp); Default: 0; Definition: exists[MulStart::ls | round(MulStart::XValue(MulStart::ls),2)=DP_x(xp)]; } StringParameter CurveSymbol { IndexDomain: xp; Definition: { if CurveMarkers(xp) then "Star" endif; } } ElementParameter CurveSymbolColor { IndexDomain: xp; Range: AllColors; Definition: { if CurveMarkers(xp) and round(x,2) = DP_x(xp) then 'red' else 'navy blue' endif; } } StringParameter FunctionFormat { Definition: FormatString("y=sin(%n*x) * sin(15*x) * (x - 2.5)",ARandomNumber); } StringParameter StatusInfo { Definition: { if SelectedPoint <> '' then FormatString("x=%n, y=%n",DP_x(SelectedPoint),DP_y(SelectedPoint)) else "" endif; } } } DeclarationSection Gantt_Chart_Declaration { Parameter GC_ClusterPoint { IndexDomain: (MulStart::cp); Definition: MulStart::cp in MulStart::ClusterPoints; } Parameter GC_ClusterStartPoint { IndexDomain: (MulStart::cp); Text: "Clusters and its radius"; Definition: MulStart::ClusterRadiusCenter(MulStart::cp) - MulStart::ClusterRadius(MulStart::cp); } Parameter GC_ClusterLength { IndexDomain: (MulStart::cp); Definition: 2 * MulStart::ClusterRadius(MulStart::cp); } StringParameter GC_ObjectiveValue { IndexDomain: (MulStart::cp); } StringParameter GC_ClusterPointDescription { IndexDomain: (MulStart::cp); Definition: "Cluster " + ord(MulStart::cp, MulStart::ClusterPoints); } } DeclarationSection Network_Declaration { Parameter NW_StartPoint_x { IndexDomain: (as); Definition: 100 * StartPoint(as); } Parameter NW_StartPoint_y { IndexDomain: (as); Definition: 0; } Parameter NW_EndPoint_x { IndexDomain: (as); Definition: 100 * EndPoint_x(as); } Parameter NW_EndPoint_y { IndexDomain: (as); Definition: 50 + 25 * EndPoint_y(as); } Parameter NW_Line { IndexDomain: (as,as2); Definition: as = as2; } ElementParameter NW_LineColor { IndexDomain: (as); Range: AllColors; Default: 'black'; } Set LinePoints { Index: lp, lp2; Definition: data{Begin, End}; } Set LablePoints { Index: lb; Definition: data { x, y }; } Parameter LP_x { IndexDomain: (lp); Definition: 500 * [ord(lp) - 1]; } Parameter LP_y { IndexDomain: (lp); Definition: 50; } Parameter LP_Line { IndexDomain: (lp,lp2); Definition: lp = lp2 - 1; } Parameter LeftBound { Definition: -10; } Parameter RightBound { Definition: 505; } StringParameter Lables { IndexDomain: lb; Definition: data { x : "Sample Points", y : "Solutions" }; } Parameter Lable_x { IndexDomain: lb; Definition: 220; } Parameter Lable_y { IndexDomain: lb; Definition: data { x : -10, y : 120 }; } } } Procedure MainInitialization { Body: { ARandomNumber:= Uniform(5,8); empty AllVariables, MulStart::ObjectiveValue, MulStart::XValue, StartPoint,AllSolves,AllSolvesCounter; } } Procedure MainExecution { Body: { AllSolvesCounter := 1; option seed = 1234; empty AllVariables, MulStart::ObjectiveValue, MulStart::XValue, StartPoint,AllSolves,AllSolvesCounter; MulStart::IterationLimit := MaxIteration; GMP1 := GMP::Instance::Generate( MP ); MulStart::DoMultiStart(GMP1, SamplePoints, SelectedSamplePoints ); } ElementParameter GMP1 { Range: AllGeneratedMathematicalPrograms; } } Procedure MainTermination { Body: { return 1; } } Module Multi_Start_Module { Prefix: MulStart; Comment: { "The basic MultiStart algorithm: ------------------------------- Input: GMP corresponding to the NLP problem, NumberOfSamplePoints, NumberOfSelectedSamplePoints. 0. Set IterationCount equal to 1. 1. Generate NumberOfSamplePoints sample points from uniform distribution. Calculate the penalized objective for all sample points and select the best NumberOfSelectedSamplePoints sample points. 2. Take one of the NumberOfSelectedSamplePoints sample points. If all sample points are processed then go to step (6). 3. For all clusters, calculate the distance between the sample point and the center of the cluster. If the distance is smaller than the radius of the cluster (i.e., the sample point belongs to the cluster) then go to step (2). 4. Solve NLP by using the sample point as its starting point to obtain a candidate local solution. For all clusters, calculate the distance between the candidate local solution and the local solution belonging to the cluster. If the distance is 0 then the candidate local solution is the same as the local solution belonging to the cluster; update the center and radius of the cluster by using the sample point and go to step (2). 5. Construct new cluster by using the mean of the sample point and the candidate local solution as its center with radius equal to half the distance between these two points. Assign the candidate local solution as the local solution belonging to the cluster. Go to step (2). 6. Increment IterationCount. If the number of iterations exceeds the IterationLimit, then go to step (7). Else go to step (1). 7. Order the local solutions and store the NumberOfBestSolutions solutions in the solution repository. The algorithm starts in fact by solving the NLP once by using the starting point (i.e., doing a normal solve). The starting point and the solution obtained are used to create, as in step (5), the first cluster." } DeclarationSection MultiStart_Control_Declarations { Comment: { "Contains parameters (\'options\') that influence the performance of the MultiStart algorithm." } Parameter NumberOfSamplePoints { Range: integer; Comment: { "The number of randomly generated points in each iteration. The sample point are randomly generated by using the intervals defined by the lower and upper bounds of the variables." } } Parameter NumberOfSelectedSamplePoints { Range: integer; Comment: { "The number of points that are selected from the set of randomly generated points in each iteration. The points that are selected are the points with the best penalty function value. If NumberOfSelectedSamplePoints equals NumberOfSamplePoints then all sample points are automatically selected." } } Parameter IterationLimit { Range: integer; InitialData: 5; Comment: "Limits the number of iterations by the MultiStart algorithm."; } Parameter NumberOfBestSolutions { Range: { {1..1000} } InitialData: 1; Comment: { "The MultiStart algorithm will return with the best \'NumberOfBestSolutions\' solutions stored at the first \'NumberOfBestSolutions\' positions in the solution repository." } } Parameter InitialPenaltyWeight { Range: nonnegative; InitialData: 1000; Comment: { "To calculate the penalty function, the violation of each constraint is multiplied by the shadow price of each constraint. This initial penalty weight acts as a lower bound on each shadow price, i.e., if the shadow price of a constraint is lower than the initial penalty weight then the initial penalty weight is used instead of the shadow price." } } Parameter ClusterRadiusMultiplier { Range: nonnegative; InitialData: 0.6; Comment: "If a new cluster is created then the radius will be multiplied by this parameter."; } Parameter ShrinkFactor { Range: [0, 1]; InitialData: 0.95; Comment: { "At the end of each iteration each cluster is multiplied by this factor to shrink the size of the cluster. Points are then less often assigned to clusters that were initially too large." } } Parameter UsePresolver { Range: binary; InitialData: 0; Comment: { "If 1, the MultiStart algorithm will start by generating the presolved GMP of the original GMP, and use that presolved GMP in the algorithm. By presolving the model, many ranges of the columns become smaller, and rows and columns might be deleted. Using this smaller GMP might speed up the MultiStart algorithm." } } Parameter UseConstraintConsensusMethod { Range: binary; InitialData: 0; Comment: { "The sample point are randomly generated by using the intervals defined by the lower and upper bounds of the variables. Such a sample point is very likely to be infeasible with respect to the constraints. The constraint consensus method (CCM) tries to find an approximately feasible point for a sample point. If 1, the MultiStart algorithm will apply CCM to each sample point generated. This might slow down the MultiStart algorithm but the chance of generating (almost) feasible sample points increases. The constraint consensus method is described in: John W. Chinneck, The Constraint Consensus Method for Finding Approximately Feasible Points in Nonlinear Programs, INFORMS Journal on Computing 16(3) (2004), pp. 255-265." } } Parameter MaximalVariableBound { Range: nonnegative; InitialData: 1000; Comment: { "If a variable has no upper bound then the upper bound used for generating sample points will be equal to the value of this parameter. If a variable has no lower bound then the value of this parameter multiplied by -1 will be used as lower bound." } } Parameter Epsilon { Definition: 1e-6; Comment: { "If the distance between two feasible solutions is smaller or equal to the value of this parameter, then the solutions are considered to be identical." } } } DeclarationSection MultiStart_Status_Declarations { Comment: "Identifiers used by the MultiStart algorithm."; ElementParameter GNLP { Range: AllGeneratedMathematicalPrograms; Comment: { "The GMP used by the MultiStart algorithm. Either the original GMP or the presolved GMP (if parameter UsePresolver equals 1)." } } ElementParameter ssNLP { Range: AllSolverSessions; Comment: "Solver session belonging to GNLP."; } Parameter OptimizationDirection { Range: { {-1..1} } Comment: "The optimization direction: 1 if maximization and -1 otherwise."; } Parameter IterationCount { Range: { {0..inf} } Property: Integer; Comment: "Iteration count for the MultiStart algorithm."; } Parameter MultiStartAlgorithmFinished { Range: binary; InitialData: 1; Comment: "Indicates whether the MultiStart algorithm has finished."; } Set Points { SubsetOf: Integers; Index: p; Comment: { "Set of integer references to points used by the MultiStart algorithm. Contains sampled points, cluster points and local solutions (as given by th sets SampledPoints, ClusterPoints and LocalSolutions respectively)." } } Set SampledPoints { SubsetOf: Points; Index: sp; Comment: { "Set of integer references to the sampled points. The first sample point is stored at position 1 and the last sample point at position NumberOfSamplePoints." } } Set ClusterPoints { SubsetOf: Points; Index: cp; Comment: { "Set of integer references to the cluster points. For each cluster two points are stored, namely the cluster centre and the local solution belonging to the cluster. These two points are stored in pairs, i.e., at position i the centre point and at position i+1 the local solution. The first cluster point is stored at position FirstClusterPoint and the last cluster point at position LastClusterPoint-1. The set LocalSolutions is the set of integer references to the local solutions." } } Set LocalSolutions { SubsetOf: Points; Index: ls; Comment: { "Set of integer references to the local solutions. For each cluster two points are stored, namely the cluster centre and the local solution belonging to the cluster. These two points are stored in pairs, i.e., at position i the centre point and at position i+1 the local solution. The first local solution is stored at position FirstClusterPoint+1 and the last local solution at position LastClusterPoint. ClusterPoints is the set of integer references to the cluster points. Local solutions are solutions returned by the NLP solver. The local solution can refer to an infeasible point." } } Parameter PenaltyWeightsPoint { Range: integer; Comment: { "The penalty function is calculated as the objective value plus the sum over all constraints of the violation multiplied by its shadow price. The shadow prices of the constraints are stored in a special point, and this parameter is the integer reference to that point." } } Parameter FirstClusterPoint { Range: integer; Comment: { "Integer reference to the first cluster point. Note that for each cluster a pair consisting of a centre point and a local solution is stored." } } Parameter LastClusterPoint { Range: integer; Comment: { "Integer reference to the last cluster point. Note that for each cluster a pair consisting of a centre point and a local solution is stored." } } Parameter PenaltyFunctionValue { IndexDomain: p; Comment: { "Parameter used to store the penalty function values of the sample points." } } Parameter SamplePointSelected { IndexDomain: p; Range: binary; Comment: { "Indicator for which sample points are selected. If a sample point is selected then the MultiStart algorithm will try to assign it to an existing cluster, and if it fails a NLP solve is started by using this sample point as the starting point. If NumberOfSelectedSamplePoints equals NumberOfSamplePoints then all sample points are automatically selected." } } Parameter ClusterRadius { IndexDomain: (cp); Range: nonnegative; Comment: "The radius of all clusters."; } Parameter ClusterRadiusInitial { IndexDomain: (cp); Range: nonnegative; Comment: "The radius of all clusters."; } Parameter ClusterRadiusCenter { IndexDomain: (cp); } Parameter NumberOfClusterPoints { IndexDomain: (cp); Range: integer; Comment: { "Number of points found by the algorithm that belong to a certain cluster. This includes the local solution corresponding to the cluster." } } Parameter ObjectiveValue { IndexDomain: (ls); Comment: { "The objective value for each local solution. Local solutions are solutions returned by the NLP solver. The local solution can refer to an infeasible point and in that case the value of the objective value is set equal to \'na\'." } } Parameter XValue { IndexDomain: (ls); } Parameter NumberOfLocalSolutions { Range: integer; Comment: { "The number of distinguished feasible local solutions found by the MultiStart algorithm. A local solution is a solution returned by the NLP solver." } } Set NLPOptimalityStatuses { SubsetOf: AllSolutionStates; Definition: { data { 'Optimal' , 'LocallyOptimal' , 'IntermediateNonOptimal' } } Comment: { "The program statusses which indicate that the solver has found a solution." } } ElementParameter ProgramStatus { Range: AllSolutionStates; Comment: { "Parameter for holding the program status corresponding to a solution." } } Parameter NumberOfSolves { Range: integer; Comment: { "The number of times the NLP solver is called to solve the NLP problem. Used to show progress information in the progress window." } } Parameter TotalSolverTime { Comment: { "The total CPU time used by the NLP solver to solve the NLP problem multiple times. Used to show progress information in the progress window." } } Parameter PresolveInfeasible { Range: binary; Comment: { "If 1, the (nonlinear) presolver has detected that the NLP problem is infeasible. In that case the MultiStart algorithm terminates immediately. This parameter is only relevant if the parameter UsePresolver is set to 1." } } } DeclarationSection Option_Declarations { StringParameter CurrentValueWarningNonOptimal { Comment: { "The current value of the option \'Warning non optimal\'. This option is temporary switched off by the MultiStart algorithm, and at the end its value is switched back to the option value just before the MultiStart algorithm was called." } } StringParameter CurrentValueWarningDerivativeEvaluation { Comment: { "The current value of the option \'Warning derivative evaluation errors\'. This option is temporary switched off by the MultiStart algorithm, and at the end its value is switched back to the option value just before the MultiStart algorithm was called." } } StringParameter CurrentValueNonlinearPresolve { Comment: { "The current value of the option \'Nonlinear presolve\'. This option is temporary switched off by the MultiStart algorithm, and at the end its value is switched back to the option value just before the MultiStart algorithm was called. Using the nonlinear presolver each time the NLP model is solved during the MultiStart algorithm is not efficient. It is much more efficient to call the presolver just once at start of the algorithm, as controlled by the parameter UsePresolver." } } Parameter OptionsChanged { Range: binary; Comment: { "If 1, the value of one of the three options was changed by the MultiStart algorithm." } } } Procedure DoMultiStart { Arguments: (MyGMP,MyNumberOfSamplePoints,MyNumberOfSelectedSamplePoints); Body: { InitializeAlgorithm( MyGMP, MyNumberOfSamplePoints, MyNumberOfSelectedSamplePoints ); if ( not PresolveInfeasible ) then DoLocalSearchFromInitialPoint; endif; while ( not MultiStartAlgorithmFinished ) do SamplePoints; AssignPointsToClusters; PrepareForNextIteration; endwhile; } Comment: "The main procedure used by the MultiStart algorithm showing the major steps."; ElementParameter MyGMP { Range: AllGeneratedMathematicalPrograms; Property: Input; } Parameter MyNumberOfSamplePoints { Range: integer; Property: Input; } Parameter MyNumberOfSelectedSamplePoints { Range: integer; Property: Input; } } Procedure InitializeAlgorithm { Arguments: (MyGMP,MyNumberOfSamplePoints,MyNumberOfSelectedSamplePoints); Body: { if ( MyNumberOfSamplePoints < MyNumberOfSelectedSamplePoints ) then halt with "Procedure DoMultiStart: 'MyNumberOfSamplePoints' smaller than 'MyNumberOfSelectedSamplePoints'."; return; endif; ChangeOptionSettings; PresolveInfeasible := 0; if ( UsePresolver ) then GNLP := GMP::Instance::CreatePresolved( MyGMP, "PresolvedModel" ); if ( GNLP = "" ) then PresolveInfeasible := 1; return AlgorithmTerminate; endif; else GNLP := MyGMP; endif; ssNLP := GMP::Instance::CreateSolverSession( GNLP ) ; NumberOfSamplePoints := MyNumberOfSamplePoints; NumberOfSelectedSamplePoints := MyNumberOfSelectedSamplePoints; IterationCount := 1; MultiStartAlgorithmFinished := 0; ! Get optimization direction. if ( GMP::Instance::GetDirection( GNLP ) = 'maximize' ) then OptimizationDirection := 1; else OptimizationDirection := -1; endif; PenaltyWeightsPoint := NumberOfSamplePoints + 1; FirstClusterPoint := NumberOfSamplePoints + 2; LastClusterPoint := FirstClusterPoint - 1; empty Points, SampledPoints, ClusterPoints; cleandependents Points, SampledPoints, ClusterPoints; empty PenaltyFunctionValue, SamplePointSelected, ClusterRadius, ObjectiveValue, XValue; while ( LoopCount <= NumberOfSamplePoints ) do Points += LoopCount; endwhile; SampledPoints := Points; NumberOfSolves := 0; NumberOfLocalSolutions := 0; TotalSolverTime := 0; ! Assign initial penalty weights to the shadow prices. GMP::Solution::UpdatePenaltyWeights( GNLP, PenaltyWeightsPoint, PenaltyWeightsPoint, InitialPenaltyWeight ); } Comment: { "Initialize the MultiStart algorithm. The model is presolved if the corresponding parameter is switched on. (Step 0 in the basic algorithm.)" } ElementParameter MyGMP { Range: AllGeneratedMathematicalPrograms; Property: Input; } Parameter MyNumberOfSamplePoints { Range: integer; Property: Input; } Parameter MyNumberOfSelectedSamplePoints { Range: integer; Property: Input; } } Procedure DoLocalSearchFromInitialPoint { Body: { CurrentPoint := 1; GMP::Solution::RetrieveFromModel( GNLP, CurrentPoint ); DoLocalSearch( CurrentPoint ); } Comment: "Call the NLP solver for the initial point as provided by the user."; Parameter CurrentPoint { Range: integer; } } Procedure SamplePoints { Body: { NumberOfPointsSelected := 0; CurrentPoint := 1; empty SelectedPoints; while ( CurrentPoint <= NumberOfSamplePoints ) do GMP::Solution::RandomlyGenerate( GNLP, CurrentPoint, MaximalVariableBound ); if ( UseConstraintConsensusMethod ) then GMP::Instance::FindApproximatelyFeasibleSolution( GNLP, CurrentPoint, CurrentPoint, IterCCM, 1000, 1e-3, 1e-3, 0.0 ); endif; PenaltyFunctionValue(CurrentPoint) := GMP::Solution::GetPenalizedObjective( GNLP, CurrentPoint, PenaltyWeightsPoint ); if ( PenaltyFunctionValue(CurrentPoint) <> na ) then SelectedPoints += CurrentPoint; NumberOfPointsSelected += 1; endif; CurrentPoint += 1; endwhile; if ( NumberOfSelectedSamplePoints <= NumberOfSamplePoints and NumberOfSelectedSamplePoints <= NumberOfPointsSelected ) then BestPoints := NBest( spts, OptimizationDirection * PenaltyFunctionValue(spts), NumberOfSelectedSamplePoints ); endif; SamplePointSelected(sp) := 0; SamplePointSelected(bp) := 1; } Comment: { "Generate NumberOfSamplePoints sample points from the uniform distribution. Calculate the penalized objective for all sample points and select the best NumberOfSelectedSamplePoints sample points. Optionally the algorithm will try to find an approximately feasible point for each sample point. (Step 1 in the basic algorithm.)" } Parameter CurrentPoint { Range: integer; } Parameter NumberOfPointsSelected { Range: integer; } Set SelectedPoints { SubsetOf: SampledPoints; Index: spts; } Set BestPoints { SubsetOf: SelectedPoints; Index: bp; OrderBy: user; } Parameter IterCCM; } Procedure AssignPointsToClusters { Body: { for ( sp | SamplePointSelected(sp) ) do ! Check whether current sample point is inside an existing cluster. AssignedToCluster := 0; for ( cp ) do Distance := GMP::Solution::GetDistance( GNLP, cp, sp ); if ( Distance <= ClusterRadius(cp) ) then if ShowAssignedPoints then AllSolvesCounter += 1; AllSolves += AllSolvesCounter; StartPoint(AllSolvesCounter) := GMP::Solution::GetColumnValue(GNLP, sp, x); EndPoint_x(AllSolvesCounter) := GMP::Solution::GetColumnValue(GNLP, cp+1, x); EndPoint_y(AllSolvesCounter) := GMP::Solution::GetObjective(GNLP, cp+1); if ShowProgress then PageRefreshAll; delay(1); endif; endif; AssignedToCluster := 1; break; endif; endfor; if ( not AssignedToCluster ) then DoLocalSearch( sp ); endif; endfor; } Comment: { "Assign selected sample points to the clusters. For all clusters, calculate the distance between the sample point and the center of the cluster. If the distance is smaller than the radius of the cluster (i.e., the sample point belongs to the cluster) then the next sample point is taken. Else solve the NLP by using the sample point as its starting point to obtain a candidate local solution. (Steps 2 and 3 in the basic algorithm.)" } Parameter AssignedToCluster { Range: binary; } Parameter Distance { Range: nonnegative; } } Procedure DoLocalSearch { Arguments: (CurrentPoint); Body: { ! Solve model with CurrentPoint as the starting point. if ShowSolvedPoints then AllSolvesCounter += 1; AllSolves += AllSolvesCounter; StartPoint(AllSolvesCounter) := GMP::Solution::GetColumnValue(GNLP, CurrentPoint, x); endif; GMP::Solution::SendToSolverSession( ssNLP, CurrentPoint ); GMP::SolverSession::Execute( ssNLP ); NumberOfSolves += 1; TotalSolverTime += GMP::SolverSession::GetCPUSecondsUsed( ssNLP ); CurrentLocalSolution := LastClusterPoint + 2; GMP::Solution::RetrieveFromSolverSession( ssNLP, CurrentLocalSolution ); if ShowSolvedPoints then EndPoint_x(AllSolvesCounter) := GMP::Solution::GetColumnValue(GNLP, CurrentLocalSolution, x); EndPoint_y(AllSolvesCounter) := GMP::Solution::GetObjective(GNLP, CurrentLocalSolution); NW_LineColor(AllSolvesCounter) := 'blue'; if ShowProgress then PageRefreshAll; delay(0.1); endif; endif; ProgramStatus := GMP::Solution::GetProgramStatus( GNLP, CurrentLocalSolution ) ; if ( ProgramStatus = 'NoSolution' or ProgramStatus = 'UnknownError' ) then GMP::Solution::Delete( GNLP, CurrentLocalSolution ); else ! Check whether new local solution is already a local solution belonging to an existing cluster. LocalSolutionFoundBefore := 0; for ( cp ) do Distance := GMP::Solution::GetDistance( GNLP, cp+1, CurrentLocalSolution ); if ( Distance <= Epsilon ) then ! CurrentPoint leads to a local solution that is already the center of an existing cluster. UpdateCluster( CurrentPoint, cp ); LocalSolutionFoundBefore := 1; endif; endfor; if ( LocalSolutionFoundBefore ) then GMP::Solution::Delete( GNLP, CurrentLocalSolution ); else ! New local solution defines a new cluster. CreateNewCluster( CurrentPoint, CurrentLocalSolution ); endif; endif; } Comment: { "Solve the NLP by using the current sample point as its starting point to obtain a candidate local solution. For all clusters, calculate the distance between the candidate local solution and the local solution belonging to the cluster. If the distance is 0 then the candidate local solution is the same as the local solution belonging to the cluster; update the center and radius of the cluster by using the sample point. (Step 4 in the basic algorithm.)" } Parameter CurrentPoint { Range: integer; Property: Input; } Parameter CurrentLocalSolution { Range: integer; } Parameter Distance { Range: nonnegative; } Parameter LocalSolutionFoundBefore { Range: binary; } } Procedure CreateNewCluster { Arguments: (CurrentPoint,CurrentLocalSolution); Body: { ! Add center of new cluster to set of ClusterPoints. ClusterCenter := LastClusterPoint + 1; Points += ClusterCenter; ClusterPoints += ClusterCenter; GMP::Solution::Copy( GNLP, CurrentPoint, ClusterCenter ); ! Construct cluster center by taking the mean of the CurrentPoint and the CurrentLocalSolution. GMP::Solution::ConstructMean( GNLP, ClusterCenter, CurrentLocalSolution, 1 ); ClusterRadius(ClusterCenter) := ClusterRadiusMultiplier * GMP::Solution::GetDistance( GNLP, CurrentPoint, CurrentLocalSolution ); ClusterRadiusInitial(ClusterCenter) := ClusterRadius(ClusterCenter); NumberOfClusterPoints(ClusterCenter) := 2; ! Add local solution belonging to new cluster to set of LocalSolutions. Points += CurrentLocalSolution; LocalSolutions += CurrentLocalSolution; ProgramStatus := GMP::Solution::GetProgramStatus( GNLP, CurrentLocalSolution ) ; if ( ProgramStatus in NLPOptimalityStatuses ) then y := GMP::Solution::GetObjective( GNLP, CurrentLocalSolution ); ObjectiveValue(CurrentLocalSolution) := y; x := GMP::Solution::GetColumnValue(GNLP,CurrentLocalSolution,x); XValue(CurrentLocalSolution) := x; ClusterRadiusCenter(ClusterCenter) := GMP::Solution::GetColumnValue(GNLP,ClusterCenter,x); NumberOfLocalSolutions += 1; else ObjectiveValue(CurrentLocalSolution) := na; endif; LastClusterPoint += 2; !PageOpen("Overview of Solutions"); !delay(1); !PageRefreshAll; ! } Comment: { "Construct new cluster by using the mean of the current sample point and the current candidate local solution as its center with radius equal to half the distance between these two points. Assign the candidate local solution as the local solution belonging to the cluster. (Step 5 in the basic algorithm.)" } Parameter CurrentPoint { Range: integer; Property: Input; } Parameter CurrentLocalSolution { Range: integer; Property: Input; } Parameter ClusterCenter { Range: integer; } Parameter Distance { Range: nonnegative; } } Procedure PrepareForNextIteration { Body: { if ( IterationCount >= IterationLimit ) then AlgorithmTerminate; else IterationCount += 1; ClusterRadius(p) *= ShrinkFactor; endif; } Comment: { "Prepare algorithm for the next iteration. If the number of iterations exceeds the IterationLimit then terminate the algorithm. (Step 6 in the basic algorithm.)" } } Procedure AlgorithmTerminate { Body: { GC_ObjectiveValue(cp) := Formatstring("%.3n",GMP::Solution::GetObjective(GNLP,cp+1)); MultiStartAlgorithmFinished := 1; ! Switch back options that were changed during the initialization phase of the algorithm. if ( OptionsChanged ) then OptionSetString( 'warning_non_optimal', CurrentValueWarningNonOptimal ); OptionSetString( 'warning_derivative_evaluation_errors', CurrentValueWarningDerivativeEvaluation ); OptionSetString( 'nonlinear_presolve', CurrentValueNonlinearPresolve ); endif; if ( PresolveInfeasible ) then NumberOfStoredSolutions := 0; else OrderLocalSolutions( NumberOfStoredSolutions ); endif; ShowFinalProgressWindow( NumberOfStoredSolutions ); if ( UsePresolver and not PresolveInfeasible ) then GMP::Instance::Delete( GNLP ); endif; } Comment: { "Terminate the MultiStart algorithm. Order the local solutions found by the algorithm. (Step 7 in the basic algorithm.)" } Parameter NumberOfStoredSolutions { Range: integer; } } Section SubRoutines { Procedure ChangeOptionSettings { Body: { ! Showing warnings for infeasible problems is unwanted. rval := OptionGetString( 'warning_non_optimal', CurrentValueWarningNonOptimal ); OptionSetString( 'warning_non_optimal', 'off' ); ! Do not show derivative evaluation errors. rval := OptionGetString( 'warning_derivative_evaluation_errors', CurrentValueWarningDerivativeEvaluation ); OptionSetString( 'warning_derivative_evaluation_errors', 'off' ); ! Switch off option for nonlinear presolve. Use control parameter 'UsePresolver' to use the ! nonlinear presolver in combination with multistart. rval := OptionGetString( 'nonlinear_presolve', CurrentValueNonlinearPresolve ); OptionSetString( 'nonlinear_presolve', 'off' ); OptionsChanged := 1; } Comment: { "Temporary switch off some options. At the end of the algorithm the options are switched back to their initial value." } Parameter rval; } Procedure UpdateCluster { Arguments: (CurrentPoint,CurrentCluster); Body: { ! Temporary copy the current cluster center to a solution buffer that is not used yet. CurrentClusterCopy := LastClusterPoint + 1; GMP::Solution::Copy( GNLP, CurrentCluster, CurrentClusterCopy ); ! Replace the cluster center by the weighted mean of the cluster center and the CurrentPoint. MeanWeight := NumberOfClusterPoints(CurrentCluster); GMP::Solution::ConstructMean( GNLP, CurrentCluster, CurrentPoint, MeanWeight ); ! Update the radius of the cluster. ClusterRadius(CurrentCluster) += GMP::Solution::GetDistance( GNLP, CurrentCluster, CurrentClusterCopy ); NumberOfClusterPoints(CurrentCluster) += 1; ! Delete temporary solution buffer used to store copy of cluster center. GMP::Solution::Delete( GNLP, CurrentClusterCopy ); } Comment: { "Update the center and radius of the current cluster by using the current sample point." } Parameter CurrentPoint { Range: integer; Property: Input; } Parameter CurrentCluster { Range: integer; Property: Input; } Parameter CurrentClusterCopy { Range: integer; } Parameter MeanWeight; } Procedure OrderLocalSolutions { Arguments: (NumberOfStoredSolutions); Body: { ! Delete sampled points. for ( sp ) do GMP::Solution::Delete( GNLP, sp ); endfor; if ( NumberOfBestSolutions < NumberOfLocalSolutions ) then NumberOfStoredSolutions := NumberOfBestSolutions; else NumberOfStoredSolutions := NumberOfLocalSolutions; endif; FeasibleLocalSolutions := { ls | ObjectiveValue(ls) <> na }; BestPoints := NBest( fls, OptimizationDirection * ObjectiveValue(fls), NumberOfStoredSolutions ); if ( NumberOfStoredSolutions <= NumberOfSamplePoints ) then CurrentLocalSolution := 1; for ( bp ) do GMP::Solution::Move( GNLP, bp, CurrentLocalSolution ); CurrentLocalSolution += 1; endfor; else CurrentLocalSolution := LastClusterPoint + 1; for ( bp ) do GMP::Solution::Move( GNLP, bp, CurrentLocalSolution ); CurrentLocalSolution += 1; endfor; CurrentLocalSolution := 1; for ( bp ) do GMP::Solution::Move( GNLP, CurrentLocalSolution+LastClusterPoint, CurrentLocalSolution ); CurrentLocalSolution += 1; endfor; endif; ! Delete solutions. SolutionSet := GMP::Solution::GetSolutionsSet( GNLP ); for ( sols | sols > NumberOfStoredSolutions ) do GMP::Solution::Delete( GNLP, sols ); endfor; ! Send best solution found to AIMMS identifiers. if ( NumberOfStoredSolutions > 0 ) then GMP::Solution::SendToModel( GNLP, 1 ); endif; } Comment: { "Order the local solutions found by the MultiStart algorithm. The solutions with the best objective values are placed at the first positions in the solution repository." } Parameter NumberOfStoredSolutions { Range: integer; Property: Output; } Set FeasibleLocalSolutions { SubsetOf: Points; Index: fls; } Set BestPoints { SubsetOf: FeasibleLocalSolutions; Index: bp; OrderBy: user; } Set SolutionSet { SubsetOf: Integers; Index: sols; } Parameter CurrentLocalSolution { Range: integer; } } Procedure ShowFinalProgressWindow { Arguments: (NumberOfStoredSolutions); Body: { if ( PresolveInfeasible ) then GMP::ProgressWindow::DisplaySolver( "Presolve" ); endif; GMP::ProgressWindow::DisplayLine( 1, "Phase", "MultiStart" ); ValueString := FormatString( "%i", IterationCount ); GMP::ProgressWindow::DisplayLine( 2, "Iterations", ValueString ); ValueString := FormatString( "%i", NumberOfSolves ); GMP::ProgressWindow::DisplayLine( 3, "Problems Solved", ValueString ); ValueString := FormatString( "%i", NumberOfLocalSolutions ); GMP::ProgressWindow::DisplayLine( 4, "Local Solutions", ValueString ); if ( NumberOfStoredSolutions = 0 ) then ValueString := FormatString( "na" ); else BestSolution := GMP::Solution::GetObjective( GNLP, 1 ); ValueString := FormatString( "%.5n", BestSolution ); endif; GMP::ProgressWindow::DisplayLine( 5, "Best Solution", ValueString ); ValueString := FormatString( "%.2n sec", TotalSolverTime / 100 ); GMP::ProgressWindow::DisplayLine( 6, "Solver Time", ValueString ); if ( PresolveInfeasible ) then GMP::ProgressWindow::DisplayProgramStatus( 'Infeasible', lineNo : 7 ) ; GMP::ProgressWindow::DisplaySolverStatus( 'NormalCompletion', lineNo : 8 ) ; elseif ( NumberOfStoredSolutions = 0 ) then GMP::ProgressWindow::DisplayProgramStatus( 'LocallyInfeasible', lineNo : 7 ) ; GMP::ProgressWindow::DisplaySolverStatus( 'NormalCompletion', lineNo : 8 ) ; else GMP::ProgressWindow::DisplayProgramStatus( 'LocallyOptimal', lineNo : 7 ) ; GMP::ProgressWindow::DisplaySolverStatus( 'NormalCompletion', lineNo : 8 ) ; endif; } Comment: { "Display the final information in the progress window at the end of the algorithm." } Parameter NumberOfStoredSolutions { Property: Input; } StringParameter ValueString; Parameter BestSolution; } } } }