ipopt-src 0.2.3+3.14.16

Redistribution of Coin-OR Ipopt as a crate
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126

// Copyright (C) 2004, 2007 International Business Machines and others.
// All Rights Reserved.
// This code is published under the Eclipse Public License.
//
// Authors:  Carl Laird, Andreas Waechter     IBM    2005-02-11

#include "IpRestoRestoPhase.hpp"
#include "IpRestoIpoptNLP.hpp"

namespace Ipopt
{
#if IPOPT_VERBOSITY > 0
static const Index dbg_verbosity = 0;
#endif

RestoRestorationPhase::RestoRestorationPhase()
{ }

RestoRestorationPhase::~RestoRestorationPhase()
{ }

bool RestoRestorationPhase::InitializeImpl(
   const OptionsList& /*options*/,
   const std::string& /*prefix*/
)
{
   return true;
}

bool RestoRestorationPhase::PerformRestoration()
{
   DBG_START_METH("RestoRestorationPhase::PerformRestoration",
                  dbg_verbosity);
   Jnlst().Printf(J_DETAILED, J_MAIN,
                  "Performing second level restoration phase for current constriant violation %8.2e\n",
                  IpCq().curr_constraint_violation());

   DBG_ASSERT(IpCq().curr_constraint_violation() > 0.);

   // Get a grip on the restoration phase NLP and obtain the pointers
   // to the original NLP data
   SmartPtr<RestoIpoptNLP> resto_ip_nlp = static_cast<RestoIpoptNLP*>(&IpNLP());
   DBG_ASSERT(dynamic_cast<RestoIpoptNLP*> (&IpNLP()));
   SmartPtr<IpoptNLP> orig_ip_nlp = static_cast<IpoptNLP*>(&resto_ip_nlp->OrigIpNLP());
   DBG_ASSERT(dynamic_cast<IpoptNLP*> (&resto_ip_nlp->OrigIpNLP()));

   // Get the current point and create a new vector for the result
   SmartPtr<const CompoundVector> Ccurr_x = static_cast<const CompoundVector*>(GetRawPtr(IpData().curr()->x()));
   DBG_ASSERT(dynamic_cast<const CompoundVector*> (GetRawPtr(IpData().curr()->x())));
   SmartPtr<const CompoundVector> Ccurr_s = static_cast<const CompoundVector*>(GetRawPtr(IpData().curr()->s()));
   DBG_ASSERT(dynamic_cast<const CompoundVector*> (GetRawPtr(IpData().curr()->s())));
   DBG_ASSERT(Ccurr_s->NComps() == 1);
   SmartPtr<Vector> new_x = IpData().curr()->x()->MakeNew();
   SmartPtr<CompoundVector> Cnew_x = static_cast<CompoundVector*>(GetRawPtr(new_x));

   // The x values remain unchanged
   SmartPtr<Vector> x = Cnew_x->GetCompNonConst(0);
   x->Copy(*Ccurr_x->GetComp(0));

   // ToDo in free mu mode - what to do here?
   Number mu = IpData().curr_mu();

   // Compute the initial values for the n and p variables for the
   // equality constraints
   Number rho = resto_ip_nlp->Rho();
   SmartPtr<Vector> nc = Cnew_x->GetCompNonConst(1);
   SmartPtr<Vector> pc = Cnew_x->GetCompNonConst(2);
   SmartPtr<const Vector> cvec = orig_ip_nlp->c(*Ccurr_x->GetComp(0));
   SmartPtr<Vector> a = nc->MakeNew();
   SmartPtr<Vector> b = nc->MakeNew();
   a->Set(mu / (2. * rho));
   a->Axpy(-0.5, *cvec);
   b->Copy(*cvec);
   b->Scal(mu / (2. * rho));
   solve_quadratic(*a, *b, *nc);
   pc->Copy(*cvec);
   pc->Axpy(1., *nc);
   DBG_PRINT_VECTOR(2, "nc", *nc);
   DBG_PRINT_VECTOR(2, "pc", *pc);

   // initial values for the n and p variables for the inequality
   // constraints
   SmartPtr<Vector> nd = Cnew_x->GetCompNonConst(3);
   SmartPtr<Vector> pd = Cnew_x->GetCompNonConst(4);
   SmartPtr<Vector> dvec = pd->MakeNew();
   dvec->Copy(*orig_ip_nlp->d(*Ccurr_x->GetComp(0)));
   dvec->Axpy(-1., *Ccurr_s->GetComp(0));
   a = nd->MakeNew();
   b = nd->MakeNew();
   a->Set(mu / (2. * rho));
   a->Axpy(-0.5, *dvec);
   b->Copy(*dvec);
   b->Scal(mu / (2. * rho));
   solve_quadratic(*a, *b, *nd);
   pd->Copy(*dvec);
   pd->Axpy(1., *nd);
   DBG_PRINT_VECTOR(2, "nd", *nd);
   DBG_PRINT_VECTOR(2, "pd", *pd);

   // Now set the trial point to the solution of the restoration phase
   // s and all multipliers remain unchanged
   SmartPtr<IteratesVector> new_trial = IpData().curr()->MakeNewContainer();
   new_trial->Set_x(*new_x);
   IpData().set_trial(new_trial);

   IpData().Append_info_string("R");

   return true;
}

void RestoRestorationPhase::solve_quadratic(
   const Vector& a,
   const Vector& b,
   Vector&       v
)
{
   v.Copy(a);
   v.ElementWiseMultiply(a);

   v.Axpy(1., b);
   v.ElementWiseSqrt();

   v.Axpy(1., a);
}

} // namespace Ipopt