1 // Copyright (C) 2005, 2006 International Business Machines and others.
2 // All Rights Reserved.
3 // This code is published under the Common Public License.
4 //
5 // $Id: IpPDPerturbationHandler.hpp 759 2006-07-07 03:07:08Z andreasw $
6 //
7 // Authors:  Carl Laird, Andreas Waechter              IBM    2005-08-04
8 
9 #ifndef __IPPDPERTURBATIONHANDLER_HPP__
10 #define __IPPDPERTURBATIONHANDLER_HPP__
11 
12 #include "IpAlgStrategy.hpp"
13 
14 namespace SimTKIpopt
15 {
16 
17   /** Class for handling the perturbation factors delta_x, delta_s,
18    *  delta_c, and delta_d in the primal dual system.  This class is
19    *  used by the PDFullSpaceSolver to handle the cases where the
20    *  primal-dual system is singular or has the wrong inertia.  The
21    *  perturbation factors are obtained based on simple heuristics,
22    *  taking into account the size of previous perturbations.
23    */
24   class PDPerturbationHandler: public AlgorithmStrategyObject
25   {
26   public:
27     /**@name Constructors/Destructors */
28     //@{
29     /** Default Constructor */
30     PDPerturbationHandler();
31     /** Default destructor */
~PDPerturbationHandler()32     virtual ~PDPerturbationHandler()
33     {}
34     //@}
35 
36     /* overloaded from AlgorithmStrategyObject */
37     virtual bool InitializeImpl(const OptionsList& options,
38                                 const std::string& prefix) override;
39 
40     /** This method must be called for each new matrix, and before any
41      *  other method for generating perturbation factors.  Usually,
42      *  the returned perturbation factors are zero, but if the system
43      *  is thought to be structurally singular, they might be
44      *  positive.  If the return value is false, no suitable
45      *  perturbation could be found. */
46     bool ConsiderNewSystem(Number& delta_x, Number& delta_s,
47                            Number& delta_c, Number& delta_d);
48 
49     /** This method returns perturbation factors for the case when the
50      *  most recent factorization resulted in a singular matrix. If
51      *  the return value is false, no suitable perturbation could be
52      *  found. */
53     bool PerturbForSingularity(Number& delta_x, Number& delta_s,
54                                Number& delta_c, Number& delta_d);
55 
56     /** This method returns perturbation factors for the case when the
57      *  most recent factorization resulted in a matrix with an
58      *  incorrect number of negative eigenvalues. If the return value
59      *  is false, no suitable perturbation could be found. */
60     bool PerturbForWrongInertia(Number& delta_x, Number& delta_s,
61                                 Number& delta_c, Number& delta_d);
62 
63     /** Just return the perturbation values that have been determined
64      *  most recently */
65     void CurrentPerturbation(Number& delta_x, Number& delta_s,
66                              Number& delta_c, Number& delta_d);
67 
68     /** Methods for IpoptType */
69     //@{
70     static void RegisterOptions(SmartPtr<RegisteredOptions> roptions);
71     //@}
72 
73   private:
74     /**@name Default Compiler Generated Methods
75      * (Hidden to avoid implicit creation/calling).
76      * These methods are not implemented and
77      * we do not want the compiler to implement
78      * them for us, so we declare them private
79      * and do not define them. This ensures that
80      * they will not be implicitly created/called. */
81     //@{
82     /** Copy Constructor */
83     PDPerturbationHandler(const PDPerturbationHandler&);
84 
85     /** Overloaded Equals Operator */
86     void operator=(const PDPerturbationHandler&);
87     //@}
88 
89     /** @name Size of the most recent non-zero perturbation. */
90     //@{
91     /** The last nonzero value for delta_x */
92     Number delta_x_last_;
93     /** The last nonzero value for delta_s */
94     Number delta_s_last_;
95     /** The last nonzero value for delta_c */
96     Number delta_c_last_;
97     /** The last nonzero value for delta_d */
98     Number delta_d_last_;
99     //@}
100 
101     /** @name Size of the most recently suggested perturbation for the
102      *  current matrix. */
103     //@{
104     /** The current value for delta_x */
105     Number delta_x_curr_;
106     /** The current value for delta_s */
107     Number delta_s_curr_;
108     /** The current value for delta_c */
109     Number delta_c_curr_;
110     /** The current value for delta_d */
111     Number delta_d_curr_;
112     //@}
113 
114     /** Flag indicating if for the given matrix the perturb for wrong
115      *  inertia method has already been called. */
116     bool get_deltas_for_wrong_inertia_called_;
117 
118     /** @name Handling structural degeneracy */
119     //@{
120     /** Type for degeneracy flags */
121     enum DegenType {
122       NOT_YET_DETERMINED,
123       NOT_DEGENERATE,
124       DEGENERATE
125     };
126 
127     /** Flag indicating whether the reduced Hessian matrix is thought
128      *  to be structurally singular. */
129     DegenType hess_degenerate_;
130 
131     /** Flag indicating whether the Jacobian of the constraints is
132      *  thought to be structurally rank-deficient. */
133     DegenType jac_degenerate_;
134 
135     /** Flag counting matrices in which degeneracy was observed in the
136      *  first successive iterations.  -1 means that there was a
137      *  non-degenerate (unperturbed) matrix at some point. */
138     Index degen_iters_;
139 
140     /** Status of current trial configuration */
141     enum TrialStatus {
142       NO_TEST,
143       TEST_DELTA_C_EQ_0_DELTA_X_EQ_0,
144       TEST_DELTA_C_GT_0_DELTA_X_EQ_0,
145       TEST_DELTA_C_EQ_0_DELTA_X_GT_0,
146       TEST_DELTA_C_GT_0_DELTA_X_GT_0
147     };
148 
149     /** Current status */
150     TrialStatus test_status_;
151     //@}
152 
153     /** @name Algorithmic parameters. */
154     //@{
155     /** Maximal perturbation for x and s. */
156     Number delta_xs_max_;
157     /** Smallest possible perturbation for x and s. */
158     Number delta_xs_min_;
159     /** Increase factor for delta_xs for first required perturbation. */
160     Number delta_xs_first_inc_fact_;
161     /** Increase factor for delta_xs for later perturbations. */
162     Number delta_xs_inc_fact_;
163     /** Decrease factor for delta_xs for later perturbations. */
164     Number delta_xs_dec_fact_;
165     /** Very first trial value for delta_xs perturbation. */
166     Number delta_xs_init_;
167     /** Size of perturbation for c and d blocks. */
168     Number delta_cd_val_;
169     /** Exponent on mu in formula for of perturbation for c and d blocks. */
170     Number delta_cd_exp_;
171     /** Flag indicating whether the new values are based on the
172      *  perturbations in the last iteration or in the more recent
173      *  iteration in which a perturbation was done. */
174     bool reset_last_;
175     /** Required number of iterations for degeneracy conclusions. */
176     Index degen_iters_max_;
177     //@}
178 
179     /** @name Auxilliary methods */
180     //@{
181     /** Internal version of PerturbForWrongInertia with the
182      *  difference, that finalize_test is not called.  Returns false
183      *  if the delta_x and delta_s parameters become too large. */
184     bool get_deltas_for_wrong_inertia(Number& delta_x, Number& delta_s,
185                                       Number& delta_c, Number& delta_d);
186 
187     /** This method is call whenever a matrix had been factorization
188      *  and is not singular.  In here, we can evaluate the outcome of
189      *  the deneracy test heuristics. */
190     void finalize_test();
191     /** Compute perturbation value for constraints */
192     Number delta_cd();
193     //@}
194 
195   };
196 
197 } // namespace Ipopt
198 
199 #endif
200