(* Modules for implementation Sequential Linear Programming *) 
(* by Zafer Gurdal *) 
(* Copyright 1995  *) 
(* Department of Engineering Science and Mechanics *) 
(* Virginia Polytechnic Institute and State University *) 
(* Blacksburg, VA 24061-0219 *) 
(* 703-231-5905 *) 
(* e-mail: gurdal@gurdal.esm.vt.edu (NeXT Mail Accepted) *)  
  
BeginPackage["SeqLinProgramming`"] 
 
SequentialLP::usage = "The following modules are available in the package 
SequentialLP         
		SLP
		SLPCenters 
		"  
  
SLP::usage = "SLP[f, {linear inequalities}, 
 {nonlinear inequalities}, {x, y,....}, {x0, y0,....} ] finds a 
 minimum of f in the domain specified by the inequalities starting
 from initial values of the variables {x0, y0,....}. The variables 
 x0, y0, ... are all assumed to be none negative.  SLP 
 employs move limits of 50% of initial design variable values, 
 with 10% reductions in move limit sizes each iteration. Move 
 limit sizes and reduction size can be specified using Options." 
    
MoveLimit::usage = "MoveLimit is an option for SLP 
 specifying the initial move limit sizes. The default setting of 
 MoveLimit is Automatic which corresponds to move limits of 50% 
 of initial design variable values. Use MoveLimit->{list} of 
 numbers, to specifiy initial move limit sizes."

ShrinkMoveLimit::usage = "ShrinkMoveLimit is an option for 
 SLP specifying the fractional reduction in move limit 
 sizes. These reductions are applied during each iteration. The 
 default setting of ShrinkMoveLimit is Automatic which corresponds 
 to 10% reduction in move limit sizes. Use ShrinkMoveLimit->some 
 fraction, to specify the fractional reduction."

SolutionHistory::usage = "SolutionHistory is an option of 
	optimizations which specifies whether the 
	entire set of cycles found will be returned in a List. 
	The default (False) is to return only the final solution."
	
PlotHistory::usage = "PlotHistory is an option 
 of optimizations which provides only the design variable values 
 for plotting purposes. Once this option is used during the solution,  
 use ProgressPlot (type ?ProgressPlot to see usage) to visualize the 
 progress of the optimization (FOR 2-D PROBLEMS ONLY).  
 The default (False) no information for plotting is returned."	
 
Epsilon::usage = "Epsilon is an option of optimizations which 
 controls the convergence criterion.  The default value of Epsilon
 is Automatic which corresponds to eps = 1.0 10^-9.  
 Use Epsilon -> eps, to change it."
	
MaxIterations::usage = "MaxIterations is an option of optimization 
 which controls the convergence criterion.  The default value of 
 MaxIterations is set to $RecursionLimit which is 256."	

Options[SLP]=
	{MoveLimit->Automatic,
	 ShrinkMoveLimit->Automatic,
	 SolutionHistory -> False, 
	 PlotHistory -> False, 
	 MaxIterations -> $RecursionLimit, 
	 Epsilon -> Automatic }

ProgressPlot::usage = "ProgressPlot[fcont,data] adds the plot 
 of solution point data onto the plot of the objective 
 function provided in fcont. For the Sequential Simplex algorithm
 use SimplexProgressPlot, instead of the ProgressPlot to visualize 
 the progress of the iterations.  
 Option AnimatePlot->True activates the generation of plots for 
 each iteration."

Options[ProgressPlot] = {AnimatePlot -> False}

Begin["`Private`"] 
 
ProgressPlot[hist_,data_,opts___?OptionQ] := 
   If[ MatrixQ[ data[[1]] ], 
       ProgressSimplex[hist, data, opts], 
       ProgressOthers[hist, data, opts] 
     ]
               
ProgressOthers[hist_,data_,opts___?OptionQ] := 
   Module[{finplot,int=data[[1]], history = hist}, 
   finplot = AnimatePlot/.{opts}/.Options[ProgressPlot];
   If[ finplot, 
       FoldList[Show[#1,Graphics[{Thickness[0.007],RGBColor[1,0,0],
                                      Line[int=Append[int,#2]]} 
                                ],
                        Graphics[{PointSize[0.02],Point[#2]}] 
                    ]&, {history}, data 
               ] ,
       Show[{history},Graphics[{Thickness[0.007],RGBColor[1,0,0],
                                      Line[data]} ], 
                       Graphics[{PointSize[0.02],Map[Point[#]&,data]}]
           ] 
      ]
         ]        

 
(* Linearize the cost function and nonlinear constraints *)
Linearize[func_,dfunc_, vars_, lpt_]:=   
  	(func /. Thread[vars->lpt]) +  
  	  Dot[(dfunc /. Thread[vars->lpt]), (vars-lpt)] //N 
       
(* Combine linearized constraints with original
   Heads after linearization *)
   
combineHead[h_,b_] :=  h[b,0] 
	
                              
SLP[f_Symbol,LnCo_List,nonLnCo_List,vars_,x0_,opts___?OptionQ] :=
            SLP[f[vars],LnCo,nonLnCo,vars,x0,opts]                              

SLP[fexp_,LnCo_List,nonLnCo_List,vars_,x0_,opts___?OptionQ] := 
  Module[ {desPt,movesize,reducerate,moves,dfexp,nlhead,nlbody,mlimits,
		dnlbody,slnH,pltH,maxIt,epep,epsCheck,lincon,linobj,
		counter=0, pltHistory={}, slnHistory={},
		epsPrecision=1.0 10^-9,test=$MaxMachineNumber}, 
(* Set move limit sizes and reduction size *)  
   {movesize,reducerate,slnH,pltH,maxIt,epep}=
      {MoveLimit,ShrinkMoveLimit,SolutionHistory,PlotHistory, 
       MaxIterations,Epsilon}/.
      {opts}/.Options[SLP];    
	If[epep === Automatic, epsCheck=epsPrecision, epsCheck=epep ];
	If[movesize===Automatic,moves=.5, moves=movesize];
	If[reducerate===Automatic,reduce=.1,reduce=reducerate];
(* Construct initial design point in the form {f(x0), {x0}} *)	
	desPt = {fexp/.Thread[vars->x0], x0};
(* Insert the appropriate quantities into appropriate history items *)	
	If[pltH, 
              pltHistory=Append[pltHistory, x0], 
              slnHistory=Append[slnHistory, desPt]     ];
(* Seperate Heads of nonlinear consts for linearization purposes *)  
	nlhead= Map[Head[#]&,nonLnCo]; 
	nlbody= Map[#[[1]] - #[[2]] &,nonLnCo]; 
(* Derivative of cost fun and nonlinear consts for linearization *)
	dfexp = Map[D[fexp,#]&,vars]; 
	dnlbody= Transpose[Map[D[nlbody,#]&,vars]];
	
(* Loop for minimizing cost function *)  
  
     While[test > epsCheck,
           prefval = desPt[[1]];  
           xi = desPt[[2]]; 
(* Linearize cost fuction and nonlinear constraints *)
   	       linobj= Linearize[fexp, dfexp, vars, xi ];
   	       lincon= Linearize[nlbody, dnlbody, vars, xi ];
(* Create move limits *) 
           If[counter===0,moves=moves,moves=moves-reduce*moves];
   	       mlimits = Map[(#<=0)&,
	          Flatten[Map[Plus[Times[#,(vars-xi)],-xi*moves]&, {1,-1}]]]; 
(* Apply ConstrainedMin to linearized system *)
   	       soln=ConstrainedMin[linobj, 
	             Join[Join[LnCo,MapThread[combineHead,{nlhead,lincon}]
	                      ], mlimits
	                 ],  vars ] ;
	       desPt={soln[[1]], vars /. soln[[2]] //N }; 
(* Store optimal obj and/or design variables into history arrays  *)
           If[pltH, 
              pltHistory=Append[pltHistory, desPt[[2]] ], 
              slnHistory=Append[slnHistory, desPt]
             ];                     
           test =  Abs[desPt[[1]] - prefval ];
           If[++counter >= maxIt, (test = 0.0; Print[ "Maximum 
	Number of Iterations ", maxIt, " Exceeded "]), ]
(* End of While *)	            
          ];
     If[pltH, pltHistory, If[slnH, slnHistory, desPt] ]             
(* End of Module *)   
	 ] 
 
End[] 
(* Protect[Linearize] *) 
(* SetAttributes[{Linearize}, Listable] *) 
 
EndPackage[]