Skip to content

Commit ff3c169

Browse files
committed
Merge branch 'feature/interior_point_solver' into feature/interior_point_example
2 parents c7b8744 + 36d0883 commit ff3c169

File tree

9 files changed

+173
-327
lines changed

9 files changed

+173
-327
lines changed

ocs2_ipm/include/ocs2_ipm/IpmPerformanceIndexComputation.h

Lines changed: 1 addition & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -87,7 +87,6 @@ PerformanceIndex computePerformanceIndex(const multiple_shooting::EventTranscrip
8787
* @param barrierParam : Barrier parameter of the interior point method
8888
* @param slackStateIneq : Slack variable of the state inequality constraints
8989
* @param slackStateInputIneq : Slack variable of the state-input inequality constraints
90-
* @param enableStateInequalityConstraints : Should be disabled at the initial node (i = 0)
9190
* @return Performance index for a single intermediate node
9291
*/
9392
PerformanceIndex computeIntermediatePerformance(OptimalControlProblem& optimalControlProblem, DynamicsDiscretizer& discretizer, scalar_t t,
@@ -134,11 +133,10 @@ PerformanceIndex computeEventPerformance(OptimalControlProblem& optimalControlPr
134133
* @param barrierParam : Barrier parameter of the interior point method
135134
* @param slackStateIneq : Slack variable of the state inequality constraints
136135
* @param slackStateInputIneq : Slack variable of the state-input inequality constraints
137-
* @param enableStateInequalityConstraints : Should be disabled at the initial node (i = 0)
138136
* @return Performance index of the interior point method
139137
*/
140138
PerformanceIndex toPerformanceIndex(const Metrics& metrics, scalar_t dt, scalar_t barrierParam, const vector_t& slackStateIneq,
141-
const vector_t& slackStateInputIneq, bool enableStateInequalityConstraints);
139+
const vector_t& slackStateInputIneq);
142140

143141
/**
144142
* Computes the PerformanceIndex of the interior point method based on a given Metrics at the event or terminal node.

ocs2_ipm/src/IpmHelpers.cpp

Lines changed: 6 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -37,7 +37,6 @@ namespace ipm {
3737
void condenseIneqConstraints(scalar_t barrierParam, const vector_t& slack, const vector_t& dual,
3838
const VectorFunctionLinearApproximation& ineqConstraint, ScalarFunctionQuadraticApproximation& lagrangian) {
3939
assert(barrierParam > 0.0);
40-
4140
const size_t nc = ineqConstraint.f.size();
4241
const size_t nu = ineqConstraint.dfdu.cols();
4342

@@ -95,8 +94,9 @@ vector_t retrieveSlackDirection(const VectorFunctionLinearApproximation& stateIn
9594

9695
vector_t retrieveDualDirection(scalar_t barrierParam, const vector_t& slack, const vector_t& dual, const vector_t& slackDirection) {
9796
assert(barrierParam > 0.0);
98-
vector_t dualDirection(slackDirection.size());
99-
dualDirection.array() = -(dual.array() * slackDirection.array() + (slack.array() * dual.array() - barrierParam)) / slack.array();
97+
vector_t dualDirection = dual.cwiseProduct(slack + slackDirection);
98+
dualDirection.array() -= barrierParam;
99+
dualDirection.array() /= -slack.array();
100100
return dualDirection;
101101
}
102102

@@ -108,14 +108,9 @@ scalar_t fractionToBoundaryStepSize(const vector_t& v, const vector_t& dv, scala
108108
return 1.0;
109109
}
110110

111-
scalar_t minFractionToBoundary = 1.0;
112-
vector_t fractionToBoundary = -marginRate * v.cwiseQuotient(dv);
113-
for (int i = 0; i < fractionToBoundary.size(); ++i) {
114-
if (fractionToBoundary[i] <= 0.0) {
115-
fractionToBoundary[i] = 1.0;
116-
}
117-
}
118-
return std::min(1.0, fractionToBoundary.minCoeff());
111+
const vector_t invFractionToBoundary = (-1.0 / marginRate) * dv.cwiseQuotient(v);
112+
const auto alpha = invFractionToBoundary.maxCoeff();
113+
return alpha > 0.0? std::min(1.0 / alpha, 1.0): 1.0;
119114
}
120115

121116
} // namespace ipm

ocs2_ipm/src/IpmPerformanceIndexComputation.cpp

Lines changed: 7 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -93,8 +93,11 @@ PerformanceIndex computeIntermediatePerformance(OptimalControlProblem& optimalCo
9393
scalar_t dt, const vector_t& x, const vector_t& x_next, const vector_t& u,
9494
scalar_t barrierParam, const vector_t& slackStateIneq, const vector_t& slackStateInputIneq,
9595
bool enableStateInequalityConstraints) {
96-
const auto metrics = multiple_shooting::computeIntermediateMetrics(optimalControlProblem, discretizer, t, dt, x, x_next, u);
97-
return toPerformanceIndex(metrics, dt, barrierParam, slackStateIneq, slackStateInputIneq, enableStateInequalityConstraints);
96+
auto metrics = multiple_shooting::computeIntermediateMetrics(optimalControlProblem, discretizer, t, dt, x, x_next, u);
97+
if (!enableStateInequalityConstraints) {
98+
metrics.stateIneqConstraint.clear();
99+
}
100+
return toPerformanceIndex(metrics, dt, barrierParam, slackStateIneq, slackStateInputIneq);
98101
}
99102

100103
PerformanceIndex computeEventPerformance(OptimalControlProblem& optimalControlProblem, scalar_t t, const vector_t& x,
@@ -110,15 +113,15 @@ PerformanceIndex computeTerminalPerformance(OptimalControlProblem& optimalContro
110113
}
111114

112115
PerformanceIndex toPerformanceIndex(const Metrics& metrics, scalar_t dt, scalar_t barrierParam, const vector_t& slackStateIneq,
113-
const vector_t& slackStateInputIneq, bool enableStateInequalityConstraints) {
116+
const vector_t& slackStateInputIneq) {
114117
PerformanceIndex performance = toPerformanceIndex(metrics, dt);
115118

116119
if (slackStateIneq.size() > 0) {
117120
performance.cost -= dt * barrierParam * slackStateIneq.array().log().sum();
118121
performance.equalityConstraintsSSE += dt * (toVector(metrics.stateIneqConstraint) - slackStateIneq).squaredNorm();
119122
}
120123

121-
if (slackStateInputIneq.size() > 0 && enableStateInequalityConstraints) {
124+
if (slackStateInputIneq.size() > 0) {
122125
performance.cost -= dt * barrierParam * slackStateInputIneq.array().log().sum();
123126
performance.equalityConstraintsSSE += dt * (toVector(metrics.stateInputIneqConstraint) - slackStateInputIneq).squaredNorm();
124127
}

ocs2_ipm/src/IpmSolver.cpp

Lines changed: 8 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -536,7 +536,6 @@ PerformanceIndex IpmSolver::setupQuadraticSubproblem(const std::vector<Annotated
536536
if (time[i].event == AnnotatedTime::Event::PreEvent) {
537537
// Event node
538538
auto result = multiple_shooting::setupEventNode(ocpDefinition, time[i].time, x[i], x[i + 1]);
539-
metrics[i] = multiple_shooting::computeMetrics(result);
540539
if (initializeSlackAndDualVariables) {
541540
slackStateIneq[i] =
542541
ipm::initializeSlackVariable(result.ineqConstraints.f, settings_.initialSlackLowerBound, settings_.initialSlackMarginRate);
@@ -545,6 +544,7 @@ PerformanceIndex IpmSolver::setupQuadraticSubproblem(const std::vector<Annotated
545544
slackStateInputIneq[i].resize(0);
546545
dualStateInputIneq[i].resize(0);
547546
}
547+
metrics[i] = multiple_shooting::computeMetrics(result);
548548
performance[workerId] += ipm::computePerformanceIndex(result, barrierParam, slackStateIneq[i]);
549549
dynamics_[i] = std::move(result.dynamics);
550550
stateInputEqConstraints_[i].resize(0, x[i].size());
@@ -567,7 +567,6 @@ PerformanceIndex IpmSolver::setupQuadraticSubproblem(const std::vector<Annotated
567567
const scalar_t ti = getIntervalStart(time[i]);
568568
const scalar_t dt = getIntervalDuration(time[i], time[i + 1]);
569569
auto result = multiple_shooting::setupIntermediateNode(ocpDefinition, sensitivityDiscretizer_, ti, dt, x[i], x[i + 1], u[i]);
570-
metrics[i] = multiple_shooting::computeMetrics(result);
571570
// Disable the state-only inequality constraints at the initial node
572571
if (i == 0) {
573572
result.stateIneqConstraints.setZero(0, x[i].size());
@@ -583,6 +582,7 @@ PerformanceIndex IpmSolver::setupQuadraticSubproblem(const std::vector<Annotated
583582
dualStateInputIneq[i] = ipm::initializeDualVariable(slackStateInputIneq[i], barrierParam, settings_.initialDualLowerBound,
584583
settings_.initialDualMarginRate);
585584
}
585+
metrics[i] = multiple_shooting::computeMetrics(result);
586586
performance[workerId] += ipm::computePerformanceIndex(result, dt, barrierParam, slackStateIneq[i], slackStateInputIneq[i]);
587587
multiple_shooting::projectTranscription(result, settings_.computeLagrangeMultipliers);
588588
dynamics_[i] = std::move(result.dynamics);
@@ -614,13 +614,13 @@ PerformanceIndex IpmSolver::setupQuadraticSubproblem(const std::vector<Annotated
614614
if (i == N) { // Only one worker will execute this
615615
const scalar_t tN = getIntervalStart(time[N]);
616616
auto result = multiple_shooting::setupTerminalNode(ocpDefinition, tN, x[N]);
617-
metrics[i] = multiple_shooting::computeMetrics(result);
618617
if (initializeSlackAndDualVariables) {
619618
slackStateIneq[N] =
620619
ipm::initializeSlackVariable(result.ineqConstraints.f, settings_.initialSlackLowerBound, settings_.initialSlackMarginRate);
621620
dualStateIneq[N] =
622621
ipm::initializeDualVariable(slackStateIneq[N], barrierParam, settings_.initialDualLowerBound, settings_.initialDualMarginRate);
623622
}
623+
metrics[i] = multiple_shooting::computeMetrics(result);
624624
performance[workerId] += ipm::computePerformanceIndex(result, barrierParam, slackStateIneq[N]);
625625
stateInputEqConstraints_[i].resize(0, x[i].size());
626626
stateIneqConstraints_[i] = std::move(result.ineqConstraints);
@@ -673,8 +673,11 @@ PerformanceIndex IpmSolver::computePerformance(const std::vector<AnnotatedTime>&
673673
const scalar_t dt = getIntervalDuration(time[i], time[i + 1]);
674674
const bool enableStateInequalityConstraints = (i > 0);
675675
metrics[i] = multiple_shooting::computeIntermediateMetrics(ocpDefinition, discretizer_, ti, dt, x[i], x[i + 1], u[i]);
676-
performance[workerId] += ipm::toPerformanceIndex(metrics[i], dt, barrierParam, slackStateIneq[i], slackStateInputIneq[i],
677-
enableStateInequalityConstraints);
676+
// Disable the state-only inequality constraints at the initial node
677+
if (i == 0) {
678+
metrics[i].stateIneqConstraint.clear();
679+
}
680+
performance[workerId] += ipm::toPerformanceIndex(metrics[i], dt, barrierParam, slackStateIneq[i], slackStateInputIneq[i]);
678681
}
679682

680683
i = timeIndex++;

ocs2_ipm/test/Exp0Test.cpp

Lines changed: 36 additions & 78 deletions
Original file line numberDiff line numberDiff line change
@@ -35,72 +35,16 @@ OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
3535

3636
#include "ocs2_ipm/IpmSolver.h"
3737

38+
#include <ocs2_core/constraint/LinearStateConstraint.h>
39+
#include <ocs2_core/constraint/LinearStateInputConstraint.h>
3840
#include <ocs2_core/initialization/DefaultInitializer.h>
3941
#include <ocs2_oc/test/EXP0.h>
4042

4143
using namespace ocs2;
4244

43-
class EXP0_StateIneqConstraints final : public StateConstraint {
44-
public:
45-
EXP0_StateIneqConstraints(const vector_t& xmin, const vector_t& xmax)
46-
: StateConstraint(ConstraintOrder::Linear), xmin_(xmin), xmax_(xmax) {}
47-
~EXP0_StateIneqConstraints() override = default;
48-
49-
EXP0_StateIneqConstraints* clone() const override { return new EXP0_StateIneqConstraints(*this); }
50-
51-
size_t getNumConstraints(scalar_t time) const override { return 4; }
52-
53-
vector_t getValue(scalar_t t, const vector_t& x, const PreComputation&) const override {
54-
vector_t e(4);
55-
e.head(2) = x - xmin_;
56-
e.tail(2) = xmax_ - x;
57-
return e;
58-
}
59-
60-
VectorFunctionLinearApproximation getLinearApproximation(scalar_t t, const vector_t& x, const PreComputation& preComp) const override {
61-
VectorFunctionLinearApproximation e;
62-
e.f = getValue(t, x, preComp);
63-
e.dfdx = matrix_t::Zero(4, x.size());
64-
e.dfdx.topLeftCorner(2, 2) = matrix_t::Identity(2, 2);
65-
e.dfdx.bottomRightCorner(2, 2) = -matrix_t::Identity(2, 2);
66-
return e;
67-
}
68-
69-
private:
70-
vector_t xmin_, xmax_;
71-
};
72-
73-
class EXP0_StateInputIneqConstraints final : public StateInputConstraint {
74-
public:
75-
EXP0_StateInputIneqConstraints(scalar_t umin, scalar_t umax) : StateInputConstraint(ConstraintOrder::Linear), umin_(umin), umax_(umax) {}
76-
~EXP0_StateInputIneqConstraints() override = default;
77-
78-
EXP0_StateInputIneqConstraints* clone() const override { return new EXP0_StateInputIneqConstraints(*this); }
79-
80-
size_t getNumConstraints(scalar_t time) const override { return 2; }
81-
82-
vector_t getValue(scalar_t t, const vector_t& x, const vector_t& u, const PreComputation&) const override {
83-
vector_t e(2);
84-
e << (u.coeff(0) - umin_), (umax_ - u.coeff(0));
85-
return e;
86-
}
87-
88-
VectorFunctionLinearApproximation getLinearApproximation(scalar_t t, const vector_t& x, const vector_t& u,
89-
const PreComputation& preComp) const override {
90-
VectorFunctionLinearApproximation e;
91-
e.f = getValue(t, x, u, preComp);
92-
e.dfdx = matrix_t::Zero(2, x.size());
93-
e.dfdu = (matrix_t(2, 1) << 1, -1).finished();
94-
return e;
95-
}
96-
97-
private:
98-
scalar_t umin_, umax_;
99-
};
100-
10145
TEST(Exp0Test, Unconstrained) {
102-
static constexpr size_t STATE_DIM = 2;
103-
static constexpr size_t INPUT_DIM = 1;
46+
constexpr size_t STATE_DIM = 2;
47+
constexpr size_t INPUT_DIM = 1;
10448

10549
// Solver settings
10650
const auto settings = []() {
@@ -141,8 +85,27 @@ TEST(Exp0Test, Unconstrained) {
14185
}
14286

14387
TEST(Exp0Test, Constrained) {
144-
static constexpr size_t STATE_DIM = 2;
145-
static constexpr size_t INPUT_DIM = 1;
88+
constexpr size_t STATE_DIM = 2;
89+
constexpr size_t INPUT_DIM = 1;
90+
91+
auto getStateBoxConstraint = [&](const vector_t& minState, const vector_t& maxState) {
92+
constexpr size_t numIneqConstraint = 2 * STATE_DIM;
93+
const vector_t e = (vector_t(numIneqConstraint) << -minState, maxState).finished();
94+
const matrix_t C =
95+
(matrix_t(numIneqConstraint, STATE_DIM) << matrix_t::Identity(STATE_DIM, STATE_DIM), -matrix_t::Identity(STATE_DIM, STATE_DIM))
96+
.finished();
97+
return std::make_unique<LinearStateConstraint>(std::move(e), std::move(C));
98+
};
99+
100+
auto getInputBoxConstraint = [&](scalar_t minInput, scalar_t maxInput) {
101+
constexpr size_t numIneqConstraint = 2 * INPUT_DIM;
102+
const vector_t e = (vector_t(numIneqConstraint) << -minInput, maxInput).finished();
103+
const matrix_t C = matrix_t::Zero(numIneqConstraint, STATE_DIM);
104+
const matrix_t D =
105+
(matrix_t(numIneqConstraint, INPUT_DIM) << matrix_t::Identity(INPUT_DIM, INPUT_DIM), -matrix_t::Identity(INPUT_DIM, INPUT_DIM))
106+
.finished();
107+
return std::make_unique<LinearStateInputConstraint>(std::move(e), std::move(C), std::move(D));
108+
};
146109

147110
// Solver settings
148111
const auto settings = []() {
@@ -173,14 +136,11 @@ TEST(Exp0Test, Constrained) {
173136
// add inequality constraints
174137
const scalar_t umin = -7.5;
175138
const scalar_t umax = 7.5;
176-
auto stateInputIneqConstraint = std::make_unique<EXP0_StateInputIneqConstraints>(umin, umax);
177-
problem.inequalityConstraintPtr->add("ubound", std::move(stateInputIneqConstraint));
139+
problem.inequalityConstraintPtr->add("ubound", getInputBoxConstraint(umin, umax));
178140
const vector_t xmin = (vector_t(2) << -7.5, -7.5).finished();
179141
const vector_t xmax = (vector_t(2) << 7.5, 7.5).finished();
180-
auto stateIneqConstraint = std::make_unique<EXP0_StateIneqConstraints>(xmin, xmax);
181-
auto finalStateIneqConstraint = std::make_unique<EXP0_StateIneqConstraints>(xmin, xmax);
182-
problem.stateInequalityConstraintPtr->add("xbound", std::move(stateIneqConstraint));
183-
problem.finalInequalityConstraintPtr->add("xbound", std::move(finalStateIneqConstraint));
142+
problem.stateInequalityConstraintPtr->add("xbound", getStateBoxConstraint(xmin, xmax));
143+
problem.finalInequalityConstraintPtr->add("xbound", getStateBoxConstraint(xmin, xmax));
184144

185145
const scalar_t startTime = 0.0;
186146
const scalar_t finalTime = 2.0;
@@ -196,18 +156,16 @@ TEST(Exp0Test, Constrained) {
196156
const auto primalSolution = solver.primalSolution(finalTime);
197157

198158
// check constraint satisfaction
199-
for (const auto& e : primalSolution.stateTrajectory_) {
200-
if (e.size() > 0) {
201-
ASSERT_TRUE(e.coeff(0) >= xmin.coeff(0));
202-
ASSERT_TRUE(e.coeff(1) >= xmin.coeff(1));
203-
ASSERT_TRUE(e.coeff(0) <= xmax.coeff(0));
204-
ASSERT_TRUE(e.coeff(1) <= xmax.coeff(1));
159+
for (const auto& x : primalSolution.stateTrajectory_) {
160+
if (x.size() > 0) {
161+
ASSERT_TRUE((x - xmin).minCoeff() >= 0);
162+
ASSERT_TRUE((xmax - x).minCoeff() >= 0);
205163
}
206164
}
207-
for (const auto& e : primalSolution.inputTrajectory_) {
208-
if (e.size() > 0) {
209-
ASSERT_TRUE(e.coeff(0) >= umin);
210-
ASSERT_TRUE(e.coeff(0) <= umax);
165+
for (const auto& u : primalSolution.inputTrajectory_) {
166+
if (u.size() > 0) {
167+
ASSERT_TRUE(u(0) - umin >= 0);
168+
ASSERT_TRUE(umax - u(0) >= 0);
211169
}
212170
}
213171
}

0 commit comments

Comments
 (0)