Skip to content

Commit 4d16558

Browse files
cvanaretcvanaret
committed
Computation of objective gradient in primal method
1 parent 662989e commit 4d16558

File tree

2 files changed

+42
-76
lines changed

2 files changed

+42
-76
lines changed

‎uno/ingredients/inequality_handling_methods/interior_point_methods/barrier_problems/PrimalExponentialBarrierProblem.cpp‎

Lines changed: 42 additions & 73 deletions
Original file line numberDiff line numberDiff line change
@@ -3,6 +3,7 @@
33

44
#include "PrimalExponentialBarrierProblem.hpp"
55
#include "ingredients/hessian_models/HessianModel.hpp"
6+
#include "linear_algebra/Indexing.hpp"
67
#include "optimization/Direction.hpp"
78
#include "optimization/Iterate.hpp"
89
#include "symbolic/UnaryNegation.hpp"
@@ -13,12 +14,8 @@
1314
namespace uno {
1415
PrimalExponentialBarrierProblem::PrimalExponentialBarrierProblem(const OptimizationProblem& problem, double barrier_parameter,
1516
const InteriorPointParameters& parameters):
16-
// no slacks: as many constraints as the number of equality constraints of the problem
17-
BarrierProblem(problem.model, problem.number_variables + 2*PrimalExponentialBarrierProblem::count_number_extra_variables(problem),
18-
0),
19-
problem(problem), number_extra_variables(PrimalExponentialBarrierProblem::count_number_extra_variables(problem)),
20-
barrier_parameter(barrier_parameter), parameters(parameters) {
21-
DEBUG << "The exponential barrier problem has " << this->number_variables << " variables\n";
17+
BarrierProblem(problem.model, problem.number_variables, 0),
18+
problem(problem), barrier_parameter(barrier_parameter), parameters(parameters) {
2219
}
2320

2421
void PrimalExponentialBarrierProblem::generate_initial_iterate(Iterate& initial_iterate) const {
@@ -45,43 +42,62 @@ namespace uno {
4542
// original objective gradient
4643
this->problem.evaluate_objective_gradient(iterate, evaluation_space, objective_gradient);
4744

48-
// contributions of inequality constraints
49-
std::vector<double> original_constraints(this->problem.number_constraints);
50-
this->problem.evaluate_constraints(iterate, original_constraints);
45+
// constraints
46+
std::vector<double> constraints(this->problem.number_constraints);
47+
this->problem.evaluate_constraints(iterate, constraints);
5148

52-
//const auto& constraints = evaluation_space.get_constraints();
53-
size_t gradient_index = this->problem.number_variables;
54-
// TODO we need to maintain two sets of constraint multipliers
49+
// multipliers
50+
std::vector<double> constraint_multipliers(this->problem.number_constraints);
5551
for (size_t constraint_index: Range(this->problem.number_constraints)) {
5652
const double lower_bound = this->problem.constraint_lower_bound(constraint_index);
5753
const double upper_bound = this->problem.constraint_upper_bound(constraint_index);
54+
constraint_multipliers[constraint_index] = 0.;
5855
if (is_finite(lower_bound)) {
59-
objective_gradient[gradient_index] = lower_bound - original_constraints[constraint_index] - this->barrier_parameter*
60-
std::log(iterate.multipliers.constraints[constraint_index]);
61-
++gradient_index;
56+
constraint_multipliers[constraint_index] += std::exp(-(constraints[constraint_index] - lower_bound)/this->barrier_parameter);
6257
}
6358
if (is_finite(upper_bound)) {
64-
objective_gradient[gradient_index] = original_constraints[constraint_index] - upper_bound - this->barrier_parameter*
65-
std::log(iterate.multipliers.constraints[constraint_index]);
66-
++gradient_index;
59+
constraint_multipliers[constraint_index] -= std::exp(-(upper_bound - constraints[constraint_index])/this->barrier_parameter);
6760
}
6861
}
62+
std::vector<double> lower_bound_multipliers(this->problem.number_variables);
63+
std::vector<double> upper_bound_multipliers(this->problem.number_variables);
6964
for (size_t variable_index: Range(this->problem.number_variables)) {
7065
const double lower_bound = this->problem.variable_lower_bound(variable_index);
7166
const double upper_bound = this->problem.variable_upper_bound(variable_index);
67+
lower_bound_multipliers[variable_index] = 0.;
68+
upper_bound_multipliers[variable_index] = 0.;
7269
if (is_finite(lower_bound)) {
73-
objective_gradient[gradient_index] = lower_bound - iterate.primals[variable_index] - this->barrier_parameter*
74-
std::log(iterate.multipliers.lower_bounds[variable_index]);
75-
++gradient_index;
70+
lower_bound_multipliers[variable_index] = std::exp(-(iterate.primals[variable_index] - lower_bound)/this->barrier_parameter);
7671
}
7772
if (is_finite(upper_bound)) {
78-
objective_gradient[gradient_index] = iterate.primals[variable_index] - upper_bound - this->barrier_parameter*
79-
std::log(iterate.multipliers.upper_bounds[variable_index]);
80-
++gradient_index;
73+
upper_bound_multipliers[variable_index] = -std::exp(-(upper_bound - iterate.primals[variable_index])/this->barrier_parameter);
8174
}
8275
}
76+
std::cout << "Constraint multipliers: "; print_vector(std::cout, constraint_multipliers);
77+
std::cout << "LB multipliers: "; print_vector(std::cout, lower_bound_multipliers);
78+
std::cout << "UB multipliers: "; print_vector(std::cout, upper_bound_multipliers);
79+
80+
// Jacobian
81+
const size_t number_jacobian_nonzeros = this->problem.number_jacobian_nonzeros();
82+
std::vector<int> jacobian_row_indices(number_jacobian_nonzeros);
83+
std::vector<int> jacobian_column_indices(number_jacobian_nonzeros);
84+
this->problem.compute_constraint_jacobian_sparsity(jacobian_row_indices.data(), jacobian_column_indices.data(),
85+
Indexing::C_indexing, MatrixOrder::COLUMN_MAJOR);
86+
std::vector<double> jacobian_values(number_jacobian_nonzeros);
87+
this->problem.evaluate_constraint_jacobian(iterate, jacobian_values.data());
88+
89+
for (size_t nonzero_index: Range(number_jacobian_nonzeros)) {
90+
const size_t constraint_index = static_cast<size_t>(jacobian_row_indices[nonzero_index]);
91+
const size_t variable_index = static_cast<size_t>(jacobian_column_indices[nonzero_index]);
92+
const double derivative = jacobian_values[nonzero_index];
93+
94+
objective_gradient[variable_index] -= derivative * constraint_multipliers[constraint_index];
95+
}
96+
97+
// bound constraints
98+
8399
std::cout << "PrimalExponentialBarrierProblem::evaluate_objective_gradient\n";
84-
for (size_t index: Range(gradient_index)) {
100+
for (size_t index: Range(this->problem.number_variables)) {
85101
std::cout << objective_gradient[index] << ' ';
86102
}
87103
std::cout << '\n';
@@ -96,7 +112,7 @@ namespace uno {
96112
}
97113

98114
size_t PrimalExponentialBarrierProblem::number_hessian_nonzeros(const HessianModel& hessian_model) const {
99-
return this->problem.number_hessian_nonzeros(hessian_model) + 2*this->number_extra_variables;
115+
return this->problem.number_hessian_nonzeros(hessian_model);
100116
}
101117

102118
void PrimalExponentialBarrierProblem::compute_constraint_jacobian_sparsity(int* /*row_indices*/, int* /*column_indices*/,
@@ -105,31 +121,7 @@ namespace uno {
105121

106122
void PrimalExponentialBarrierProblem::compute_hessian_sparsity(const HessianModel& hessian_model, int* row_indices,
107123
int* column_indices, int solver_indexing) const {
108-
// original Lagrangian Hessian
109124
this->problem.compute_hessian_sparsity(hessian_model, row_indices, column_indices, solver_indexing);
110-
111-
// two copies of the Jacobian
112-
const size_t number_hessian_nonzeros = this->problem.number_hessian_nonzeros(hessian_model);
113-
// sparsity of the Jacobian
114-
this->problem.compute_constraint_jacobian_sparsity(row_indices + number_hessian_nonzeros, column_indices +
115-
number_hessian_nonzeros, solver_indexing, MatrixOrder::COLUMN_MAJOR); // TODO
116-
const size_t number_jacobian_nonzeros = this->problem.number_jacobian_nonzeros();
117-
// copy it a second time
118-
for (size_t index: Range(number_jacobian_nonzeros)) {
119-
row_indices[number_hessian_nonzeros + number_jacobian_nonzeros + index] = row_indices[number_hessian_nonzeros + index];
120-
column_indices[number_hessian_nonzeros + number_jacobian_nonzeros + index] = column_indices[number_hessian_nonzeros + index];
121-
}
122-
// shift the row indices
123-
for (size_t index: Range(number_jacobian_nonzeros)) {
124-
row_indices[number_hessian_nonzeros + index] += static_cast<int>(this->problem.number_variables);
125-
row_indices[number_hessian_nonzeros + number_jacobian_nonzeros + index] += static_cast<int>(this->problem.number_variables +
126-
this->problem.number_constraints);
127-
}
128-
// diagonal contributions
129-
row_indices[11] = 7;
130-
column_indices[11] = 1;
131-
row_indices[12] = 8;
132-
column_indices[12] = 1;
133125
}
134126

135127
void PrimalExponentialBarrierProblem::evaluate_constraint_jacobian(Iterate& /*iterate*/, double* /*jacobian_values*/) const {
@@ -252,27 +244,4 @@ namespace uno {
252244
}};
253245
return norm_inf(shifted_bound_complementarity); // TODO use a generic norm
254246
}
255-
256-
size_t PrimalExponentialBarrierProblem::count_number_extra_variables(const OptimizationProblem& problem) {
257-
size_t number_variables = 0;
258-
259-
// count the number of variables associated to constraint bounds
260-
for (size_t constraint_index: Range(problem.number_constraints)) {
261-
if (is_finite(problem.constraint_lower_bound(constraint_index))) {
262-
number_variables++;
263-
}
264-
if (is_finite(problem.constraint_upper_bound(constraint_index))) {
265-
number_variables++;
266-
}
267-
}
268-
for (size_t variable_index: Range(problem.number_variables)) {
269-
if (is_finite(problem.variable_lower_bound(variable_index))) {
270-
number_variables++;
271-
}
272-
if (is_finite(problem.variable_upper_bound(variable_index))) {
273-
number_variables++;
274-
}
275-
}
276-
return number_variables;
277-
}
278247
} // namespace

‎uno/ingredients/inequality_handling_methods/interior_point_methods/barrier_problems/PrimalExponentialBarrierProblem.hpp‎

Lines changed: 0 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -70,13 +70,10 @@ namespace uno {
7070

7171
protected:
7272
const OptimizationProblem& problem;
73-
const size_t number_extra_variables;
7473
const double barrier_parameter;
7574
const InteriorPointParameters& parameters;
7675
const Vector<size_t> fixed_variables{};
7776
const ForwardRange empty_set{0};
78-
79-
[[nodiscard]] static size_t count_number_extra_variables(const OptimizationProblem& problem);
8077
};
8178
} // namespace
8279

0 commit comments

Comments
 (0)