FuncPieceError not working as expected
AnsweredDear Gurobi-Team, dear community,
I have trouble understanding how FuncPieceError works. I have compiled a minimal working of trying to piecewise linearly approximate the function
y = f(x) = x^2/(2*d)
using the general function constraints provided by Gurobi. I always want to overestimate the function value. I am ok with a somewhat large error which I try to specifically provide by setting the FuncPieceError attribute, in the example to an absolute error of 20. I also included bounds to both x and y variables in the model. I am using the C++ API. I believe to also have set FuncPieces correctly according to https://www.gurobi.com/documentation/10.0/refman/constraints.html#subsubsection:GenConstrFunction and https://www.gurobi.com/documentation/10.0/refman/funcpieceerror2.html#parameter:FuncPieceError, i.e.,
FuncPieceRatio=1 FuncPieces=-1 FuncPieceError=20
Now I am getting warning messages of the type "max general constraint violation (xxx) exceeds tolerance", this is expected since I set a large tolerance. However the specified constraint violation is sometimes larger than the specified error, e.g., 27, which is significantly larger than 20 I feel (more details see below). If I check the value of y to the actual function value, the error is eaven larger. Am I misunderstanding something? How can I specify that Gurobi should choose a piecewise representation (in the interval specified by the variable bounds) that has a certain maximal absolute error? Why does Gurobi behave as it does (see below)?
-------------------------------------
Finally, let me provide with what I have done:
I used the following C++ code:
#include <iostream>
#include "gurobi_c++.h"
#include <vector>
int main() {
std::vector<double> x_vals;
std::vector<double> y_vals;
std::vector<double> constraint_violations;
std::vector<double> real_vals;
std::vector<double> gen_constr_approximation_errors;
const double max_x = 20;
for (double target_x = 0; target_x <= max_x; target_x += 1) {
GRBEnv env = new GRBEnv(true);
env.start();
GRBModel model = GRBModel(env);
const double d = 1;
const double max_y = max_x * max_x / (2 * d);
GRBVar x = model.addVar(0, max_x, 0, GRB_CONTINUOUS, "x");
GRBVar y = model.addVar(0, max_y, 0, GRB_CONTINUOUS, "y");
GRBVar z = model.addVar(0, GRB_INFINITY, 0, GRB_CONTINUOUS, "z");
GRBLinExpr obj = z;
model.setObjective(obj, GRB_MINIMIZE);
const double p[3]{1 / (2 * d), 0, 0};
model.addGenConstrPoly(x, y, 3, p, "gen_const",
"FuncPieceRatio=1 FuncPieces=-1 FuncPieceError=20");
model.addConstr(z, GRB_GREATER_EQUAL, target_x - x, "z_constr_1");
model.addConstr(z, GRB_GREATER_EQUAL, x - target_x, "z_constr_2");
model.optimize();
if (target_x == 4) {
model.write("model.lp");
model.write("model.mps");
model.write("model.sol");
}
x_vals.emplace_back(x.get(GRB_DoubleAttr_X));
y_vals.emplace_back(y.get(GRB_DoubleAttr_X));
constraint_violations.emplace_back(model.get(GRB_DoubleAttr_ConstrVio));
real_vals.emplace_back(x_vals.back() * x_vals.back() / (2 * d));
gen_constr_approximation_errors.emplace_back(std::abs(y_vals.back() - real_vals.back()));
}
for (int i = 0; i < x_vals.size(); ++i) {
std::cout << "x=" << x_vals[i] << ", y=" << y_vals[i] << "(real: " << real_vals[i] << ", error: " << gen_constr_approximation_errors[i] << "), constraint violation: " << constraint_violations[i] << std::endl;
}
}
By this, we can conclude the following approximation used by Gurobi:
x=0, y=0(real: 0, error: 0), constraint violation: 0
x=1, y=10(real: 0.5, error: 9.5), constraint violation: 8.07587
x=2, y=20(real: 2, error: 18), constraint violation: 18
x=3, y=30(real: 4.5, error: 25.5), constraint violation: 25.5
x=4, y=40(real: 8, error: 32), constraint violation: 27.514
x=5, y=50(real: 12.5, error: 37.5), constraint violation: 25.5977
x=6, y=60(real: 18, error: 42), constraint violation: 23.086
x=7, y=70(real: 24.5, error: 45.5), constraint violation: 20.3986
x=8, y=80(real: 32, error: 48), constraint violation: 17.7254
x=9, y=90(real: 40.5, error: 49.5), constraint violation: 15.1615
x=10, y=100(real: 50, error: 50), constraint violation: 12.7583
x=11, y=110(real: 60.5, error: 49.5), constraint violation: 10.545
x=12, y=120(real: 72, error: 48), constraint violation: 8.53924
x=13, y=130(real: 84.5, error: 45.5), constraint violation: 6.75118
x=14, y=140(real: 98, error: 42), constraint violation: 5.18622
x=15, y=150(real: 112.5, error: 37.5), constraint violation: 3.84557
x=16, y=160(real: 128, error: 32), constraint violation: 2.7255
x=17, y=170(real: 144.5, error: 25.5), constraint violation: 1.81494
x=18, y=180(real: 162, error: 18), constraint violation: 1.09033
x=19, y=190(real: 180.5, error: 9.5), constraint violation: 0.508171
x=20, y=200(real: 200, error: 0), constraint violation: 0
which is for some values significantly larger than the desired error of maximum 20. Also the constraint violation specified by Gurobi and the value of f(x)-y calculated by myself does not always match, sometimes it significantly varies.
Exemplarily for x=4 I provide model details:
model.lp:
\ LP format - for model browsing. Use MPS format to capture full model detail.
Minimize
0 y + z
Subject To
z_constr_1: x + z >= 4
z_constr_2: - x + z >= -4
Bounds
x <= 20
y <= 200
General Constraints
gen_const: ( FuncPieces=-1 FuncPieceError=20 FuncPieceRatio=1 ) y = POLY (
0.5 x ^ 2 )
End
model.mps:
NAME (null)
ROWS
N OBJ
G z_constr_1
G z_constr_2
COLUMNS
x z_constr_1 1
x z_constr_2 -1
y OBJ 0
z OBJ 1
z z_constr_1 1
z z_constr_2 1
RHS
RHS1 z_constr_1 4
RHS1 z_constr_2 -4
BOUNDS
UP BND1 x 20
UP BND1 y 200
GENCONS
POLY gen_const
Options -1 0.01 20 1
x y
2 0.5
ENDATA
model.sol:
# Objective value = 0
x 4
y 40
z 0
The first lines of the Gurobi Output (including operating system, Gurobi version etc.) are
Academic license - for non-commercial use only - expires 2023-12-02
Gurobi Optimizer version 10.0.1 build v10.0.1rc0 (win64)
CPU model: 12th Gen Intel(R) Core(TM) i7-1260P, instruction set [SSE2|AVX|AVX2]
Thread count: 12 physical cores, 16 logical processors, using up to 16 threads
--------------------
Thanks in advance for any help on this problem / misunderstanding. Best,
Stefan
-
-
Gurobi Staff
Hi Stefan,
You are allowing for an approximation error of at most 20 by setting FuncPieceError = 20. This parameter should be set to something like FuncPieceError=1e-3 or a bit lower. Note that Gurobi's piecewise-linear approximation is performed in presolve. Thus, the approximation will not be recomputed during the B&B algorithm meaning that if the approximation introduces a big error, it will stay during the whole optimization process. This often leads to warnings about violation of general constraints because they only approximate the original function with a fixed error ratio.
For your function of interest, you don't have to use general constraints. It is a quadratic function, thus you can directly model it in Gurobi instead of using piecewise-linear approximation. In C++ you would use a GRBQuadExpr object and the addQConstr method. For an example of a quadratic model, please refer to bilinear_c++.cpp. For more details of how Gurobi handles quadratic terms, please refer to our NonConvex Webinar and Tech Talk.
Best regards,
Jaromił1 -
-
Hi Jarmoil,
thank you for your qick response and pointint to the addQConstr method. I have a few follow up questions though, because I have still trouble of understanding the exact meaning of FuncPieceError / it behaves differently than I would have expected from reading the docs.
I was trying to use piecewise linear functions, because I have a significant amount of binary variables in the full model. I believe(d) that MILPS are significantly easier to solve compared to mixed integer programs with linear objective and quadratic / nonlinear equality constraints (MIQCP, right?) and global optimality usually not guaranteed for the latter. I now read that Gurobi uses spatial branching, so Gurobi can solve to optimality similarly to classic MILPS, right? Do you have any feeling on if spatial branching or defining PWLs (hence classical MILP) is usually faster?
While you say that FuncPieceError should be something in the order of 1e-3 (though my modeling detail really doesn't care if the error is 1e-3 or around say 20 in practice so that error is really ok for me), I still have trouble understanding what the value does. I expected an approximation error of AT MOST 20 (and as a user I would expect this independently of if the chosen error is reasonable or not, also given that by hand I can calculate the number of pieces needed for this, detail see below).
However to me it seems that Gurobi sometimes allowes a larger error (more details in my little "experiment" below). Can you explain why that is? Or did you mean that this is because the absolute error is ensured in some altered presolved model, which might differ from the original model? But that would also not really make sense, because I (as a user) would expect my parameters hold for the constraint that I have specified.
Using WolframAlpha I calculated the number of expected pieces by hand, more precisely the maximal error is <= ABS_ERROR if the number of pieces N >= MAX_X/(2*std::sqrt(2)*std::sqrt(D)*std::sqrt(ABS_ERROR)). Hence, I would expect Gurobi to use
N = std::ceil(MAX_X/(2*std::sqrt(2)*std::sqrt(D)*std::sqrt(ABS_ERROR)));
pieces. In the code attached (I put the code at the very end of this message so that the following output makes sense). This seems to always produce the desired absolute approximation error, however Gurobi seems to only sometimes and not always use this.
Gurobi seems to use exactly this expected approximation for MAX_X=20, D=1.0, ABS_ERROR=5, i.e.,
PWL Gurobi and PWL floor are NOT equal
PWL Gurobi and PWL ceil are equal
PWL Gurobi and PWL round are NOT equalHowever for MAX_X=20, D=1.0, ABS_ERROR=10 Gurobi seems to floor the number and hence obtain an error larger than the specified 10 for some value (e.g. for x = 5 the error is 12.5 > 10), i.e.,
x = 0, f(x) = 0
PWL Gurobi: 0, error = 0
PWL ceil: 0, error = 0
PWL floor: 0, error = 0
PWL round: 0, error = 0
x = 1, f(x) = 0.5
PWL Gurobi: 5, error = 4.5
PWL ceil: 3, error = 2.5
PWL floor: 5, error = 4.5
PWL round: 5, error = 4.5
x = 2, f(x) = 2
PWL Gurobi: 10, error = 8
PWL ceil: 6, error = 4
PWL floor: 10, error = 8
PWL round: 10, error = 8
x = 3, f(x) = 4.5
PWL Gurobi: 15, error = 10.5
PWL ceil: 9, error = 4.5
PWL floor: 15, error = 10.5
PWL round: 15, error = 10.5
x = 4, f(x) = 8
PWL Gurobi: 20, error = 12
PWL ceil: 12, error = 4
PWL floor: 20, error = 12
PWL round: 20, error = 12
x = 5, f(x) = 12.5
PWL Gurobi: 25, error = 12.5
PWL ceil: 15, error = 2.5
PWL floor: 25, error = 12.5
PWL round: 25, error = 12.5
x = 6, f(x) = 18
PWL Gurobi: 30, error = 12
PWL ceil: 18, error = 0
PWL floor: 30, error = 12
PWL round: 30, error = 12
x = 7, f(x) = 24.5
PWL Gurobi: 35, error = 10.5
PWL ceil: 27.5, error = 3
PWL floor: 35, error = 10.5
PWL round: 35, error = 10.5
x = 8, f(x) = 32
PWL Gurobi: 40, error = 8
PWL ceil: 37, error = 5
PWL floor: 40, error = 8
PWL round: 40, error = 8
x = 9, f(x) = 40.5
PWL Gurobi: 45, error = 4.5
PWL ceil: 46.5, error = 6
PWL floor: 45, error = 4.5
PWL round: 45, error = 4.5
x = 10, f(x) = 50
PWL Gurobi: 50, error = 0
PWL ceil: 56, error = 6
PWL floor: 50, error = 0
PWL round: 50, error = 0
x = 11, f(x) = 60.5
PWL Gurobi: 65, error = 4.5
PWL ceil: 65.5, error = 5
PWL floor: 65, error = 4.5
PWL round: 65, error = 4.5
x = 12, f(x) = 72
PWL Gurobi: 80, error = 8
PWL ceil: 75, error = 3
PWL floor: 80, error = 8
PWL round: 80, error = 8
x = 13, f(x) = 84.5
PWL Gurobi: 95, error = 10.5
PWL ceil: 84.5, error = 0
PWL floor: 95, error = 10.5
PWL round: 95, error = 10.5
x = 14, f(x) = 98
PWL Gurobi: 110, error = 12
PWL ceil: 101, error = 3
PWL floor: 110, error = 12
PWL round: 110, error = 12
x = 15, f(x) = 112.5
PWL Gurobi: 125, error = 12.5
PWL ceil: 117.5, error = 5
PWL floor: 125, error = 12.5
PWL round: 125, error = 12.5
x = 16, f(x) = 128
PWL Gurobi: 140, error = 12
PWL ceil: 134, error = 6
PWL floor: 140, error = 12
PWL round: 140, error = 12
x = 17, f(x) = 144.5
PWL Gurobi: 155, error = 10.5
PWL ceil: 150.5, error = 6
PWL floor: 155, error = 10.5
PWL round: 155, error = 10.5
x = 18, f(x) = 162
PWL Gurobi: 170, error = 8
PWL ceil: 167, error = 5
PWL floor: 170, error = 8
PWL round: 170, error = 8
x = 19, f(x) = 180.5
PWL Gurobi: 185, error = 4.5
PWL ceil: 183.5, error = 3
PWL floor: 185, error = 4.5
PWL round: 185, error = 4.5
x = 20, f(x) = 200
PWL Gurobi: 200, error = 0
PWL ceil: 200, error = 0
PWL floor: 200, error = 0
PWL round: 200, error = 0
PWL Gurobi and PWL floor are equal
PWL Gurobi and PWL ceil are NOT equal
PWL Gurobi and PWL round are equalThank you in advance. Best regards,
Stefan
------------------------
main.cpp used:
#include <iostream>
#include "gurobi_c++.h"
#include <vector>
#include <memory>
#define MAX_X 20
#define D 1.0
#define ABS_ERROR 10
void add_standard_vars_and_constraints(GRBModel& model, const int& target_x) {
const double max_y = MAX_X * MAX_X / (2 * D);
GRBVar x = model.addVar(0, MAX_X, 0, GRB_CONTINUOUS, "x");
GRBVar y = model.addVar(0, max_y, 0, GRB_CONTINUOUS, "y");
GRBVar z = model.addVar(0, GRB_INFINITY, 0, GRB_CONTINUOUS, "z");
GRBLinExpr obj = z;
model.setObjective(obj, GRB_MINIMIZE);
model.addConstr(z, GRB_GREATER_EQUAL, target_x - x, "z_constr_1");
model.addConstr(z, GRB_GREATER_EQUAL, x - target_x, "z_constr_2");
model.update();
}
void add_pwl_constraint(GRBModel& model, const int& N) {
auto x = model.getVarByName("x");
auto y = model.getVarByName("y");
std::unique_ptr<double[]> xpts(new double[N+1]);
std::unique_ptr<double[]> ypts(new double[N+1]);
for (int i = 0; i <= N; ++i) {
xpts[i] = i * MAX_X / N;
ypts[i] = xpts[i] * xpts[i] / (2 * D);
}
model.addGenConstrPWL(x, y, N+1, xpts.get(), ypts.get(), "gen_const");
model.update();
}
int main() {
std::vector<int> x_vals;
std::vector<double> y_vals_pwl_floor;
std::vector<double> y_vals_pwl_ceil;
std::vector<double> y_vals_pwl_round;
std::vector<double> y_vals_pwl_gurobi;
std::vector<double> real_vals;
std::vector<double> gen_constr_approximation_errors_pwl_floor;
std::vector<double> gen_constr_approximation_errors_pwl_ceil;
std::vector<double> gen_constr_approximation_errors_pwl_round;
std::vector<double> gen_constr_approximation_errors_pwl_gurobi;
for (int target_x = 0; target_x <= MAX_X; target_x += 1) {
x_vals.emplace_back(target_x);
real_vals.emplace_back(x_vals.back() * x_vals.back() / (2 * D));
}
for (int target_x = 0; target_x <= MAX_X; target_x += 1) {
GRBEnv env = new GRBEnv(true);
env.start();
GRBModel model = GRBModel(env);
add_standard_vars_and_constraints(model, target_x);
auto x = model.getVarByName("x");
auto y = model.getVarByName("y");
const double p[3]{1 / (2 * D), 0, 0};
model.addGenConstrPoly(x, y, 3, p, "gen_const",
"FuncPieceRatio=1 FuncPieces=-1 FuncPieceError=" + std::to_string(ABS_ERROR));
model.update();
model.optimize();
y_vals_pwl_gurobi.emplace_back(y.get(GRB_DoubleAttr_X));
gen_constr_approximation_errors_pwl_gurobi.emplace_back(std::abs(y_vals_pwl_gurobi.back() - real_vals[target_x]));
}
for (int target_x = 0; target_x <= MAX_X; target_x += 1) {
GRBEnv env = new GRBEnv(true);
env.start();
GRBModel model = GRBModel(env);
add_standard_vars_and_constraints(model, target_x);
const int N = std::ceil(MAX_X/(2*std::sqrt(2)*std::sqrt(D)*std::sqrt(ABS_ERROR)));
add_pwl_constraint(model, N);
model.optimize();
auto y = model.getVarByName("y");
y_vals_pwl_ceil.emplace_back(y.get(GRB_DoubleAttr_X));
gen_constr_approximation_errors_pwl_ceil.emplace_back(std::abs(y_vals_pwl_ceil.back() - real_vals[target_x]));
}
for (int target_x = 0; target_x <= MAX_X; target_x += 1) {
GRBEnv env = new GRBEnv(true);
env.start();
GRBModel model = GRBModel(env);
add_standard_vars_and_constraints(model, target_x);
const int N = std::floor(MAX_X/(2*std::sqrt(2)*std::sqrt(D)*std::sqrt(ABS_ERROR)));
add_pwl_constraint(model, N);
model.optimize();
auto y = model.getVarByName("y");
y_vals_pwl_floor.emplace_back(y.get(GRB_DoubleAttr_X));
gen_constr_approximation_errors_pwl_floor.emplace_back(std::abs(y_vals_pwl_floor.back() - real_vals[target_x]));
}
for (int target_x = 0; target_x <= MAX_X; target_x += 1) {
GRBEnv env = new GRBEnv(true);
env.start();
GRBModel model = GRBModel(env);
add_standard_vars_and_constraints(model, target_x);
const int N = std::round(MAX_X/(2*std::sqrt(2)*std::sqrt(D)*std::sqrt(ABS_ERROR)));
add_pwl_constraint(model, N);
model.optimize();
auto y = model.getVarByName("y");
y_vals_pwl_round.emplace_back(y.get(GRB_DoubleAttr_X));
gen_constr_approximation_errors_pwl_round.emplace_back(std::abs(y_vals_pwl_round.back() - real_vals[target_x]));
}
bool equal_to_floor = true;
bool equal_to_ceil = true;
bool equal_to_round = true;
for (int i = 0; i < x_vals.size(); ++i) {
std::cout << "x = " << x_vals[i] << ", f(x) = " << real_vals[i] << std::endl;
std::cout << " PWL Gurobi: " << y_vals_pwl_gurobi[i] << ", error = " << gen_constr_approximation_errors_pwl_gurobi[i] << std::endl;
std::cout << " PWL ceil: " << y_vals_pwl_ceil[i] << ", error = " << gen_constr_approximation_errors_pwl_ceil[i] << std::endl;
std::cout << " PWL floor: " << y_vals_pwl_floor[i] << ", error = " << gen_constr_approximation_errors_pwl_floor[i] << std::endl;
std::cout << " PWL round: " << y_vals_pwl_round[i] << ", error = " << gen_constr_approximation_errors_pwl_round[i] << std::endl;
if (y_vals_pwl_gurobi[i] != y_vals_pwl_floor[i]) {
equal_to_floor = false;
}
if (y_vals_pwl_gurobi[i] != y_vals_pwl_ceil[i]) {
equal_to_ceil = false;
}
if (y_vals_pwl_gurobi[i] != y_vals_pwl_round[i]) {
equal_to_round = false;
}
}
if (equal_to_floor) {
std::cout << "PWL Gurobi and PWL floor are equal" << std::endl;
} else {
std::cout << "PWL Gurobi and PWL floor are NOT equal" << std::endl;
}
if (equal_to_ceil) {
std::cout << "PWL Gurobi and PWL ceil are equal" << std::endl;
} else {
std::cout << "PWL Gurobi and PWL ceil are NOT equal" << std::endl;
}
if (equal_to_round) {
std::cout << "PWL Gurobi and PWL round are equal" << std::endl;
} else {
std::cout << "PWL Gurobi and PWL round are NOT equal" << std::endl;
}
}0 -
-
-
Hi Stefan,
I believe(d) that MILPS are significantly easier to solve compared to mixed integer programs with linear objective and quadratic / nonlinear equality constraints (MIQCP, right?) and global optimality usually not guaranteed for the latter.
The algorithms and heuristics are more advanced for MILPs than for quadratic / nonlinear programs. Thus, it is most often better to solve a MILP instead of a nonlinear problem. However as usual, there are exceptions.
It is possible to guarantee global optimality for MINLPs through, e.g., the spatial B&B algorithm. If you use a local solution method then no guarantee of global optimality can be given.
I now read that Gurobi uses spatial branching, so Gurobi can solve to optimality similarly to classic MILPS, right?
Yes, see our Nonconvex Webinar and Tech Talk for more details.
Do you have any feeling on if spatial branching or defining PWLs (hence classical MILP) is usually faster?
Usually, if a suitable algorithm is available, I would always first try to solve the nonlinear model instead of directly using an approximation. A disadvantage of the piecewise-linear approximation is that often a lot of binary variables are needed to properly approximate the original nonlinear function. This bloats up the model size a lot and makes the MILP way harder to solve. Moreover, one always introduces an approximation error, which, depending on the application, may or may not be acceptable.
However to me it seems that Gurobi sometimes allowes a larger error (more details in my little "experiment" below). Can you explain why that is? Or did you mean that this is because the absolute error is ensured in some altered presolved model, which might differ from the original model? But that would also not really make sense, because I (as a user) would expect my parameters hold for the constraint that I have specified.
I am not sure I understand the issue here. Could you please share a model in LP or MPS format for which you see a constraint violation larger than 20 despite you setting FuncPieceError=20 and FuncPieces=-1. Note that uploading files in the Community Forum is not possible but we discuss an alternative in Posting to the Community Forum.
However for MAX_X=20, D=1.0, ABS_ERROR=10 Gurobi seems to floor the number and hence obtain an error larger than the specified 10 for some value (e.g. for x = 5 the error is 12.5 > 10), i.e.,
I don't understand how you got these values. In particular how you computed the Gurobi error at a given point. Could you please elaborate?
Edit: I think I now understand how you compute the error values. Could you still please write out the particular LP/MPS file where the error is violated such that I could have a look?
Note that Gurobi does not necessarily use a uniform grid to approximate a given function.
Best regards,
Jaromił0 -
-
Hi Jaromil,
thank you for you answer, I am happy to provide an mps file. Further down in my post, I will also try to explain my thoughts a bit more (hope that helps and clarifies what I mean).
In Short:
The mps file provided here: https://gigamove.rwth-aachen.de/de/download/e1bf82759cf5607bdf463c5f44b49d85 (passphrase: almighty-jujitsu-deluxe) can be downloaded until June 29 has a pwl approximation of y=0.5x^2. Optimal solution is x=4 and y=40 (and z=0). But the "real" value of the function at x=4 is 8. Hence, the pwl absolute error at x=4 is 32. Moreover gurobi prints a constraint violation of around 27.
I can also post the content of the mps file here in plain text. It is:NAME (null)
ROWS
N OBJ
G z_constr_1
G z_constr_2
COLUMNS
x z_constr_1 1
x z_constr_2 -1
y OBJ 0
z OBJ 1
z z_constr_1 1
z z_constr_2 1
RHS
RHS1 z_constr_1 4
RHS1 z_constr_2 -4
BOUNDS
UP BND1 x 20
UP BND1 y 200
GENCONS
POLY gen_const
Options -1 0.01 20 -1
x y
2 0.5
ENDATASome details on my reasoning:
Let me share my thoughts on the minimal example above. I want to approximate the function \(yf(x)=0.5x^2\) using a picewise linear approximation. In order to access the approximated value at a given x value, I want to make sure that Gurobi picks a certain target x (but I still want to remain a feasible region in a certain range). So lets say, we want to find f(4) but want gurobi to approximate the function on the interval [0,20]. Then I modelmin abs(x-4)
0 <= x <= 20
y = 0.5x^2 (using pwl, FuncPieces=-1 FuncPieceError=20)i.e., with an additional variable
min z
z >= x-4
z >= 4-x
y = 0.5x^2 (pwl)
0 <= x <= 20Gurobi returns as an optimal solution x=4 and y=40. However f(4)=8, hence the difference is 32, which is larger than 20. Gurobi specifies the constraint violaton as 27.514 (I am not sure what this value means in this context though, since I feel like the actual error is 32, but both values are significantly above 20). I have saved the mps file (as written above). I provide the cpp code at the bottom of this post.
Pleas note that I am getting the same results (x=4, y=40) when running gurobi from the command line as follows:gurobi_cl ResultFile=model_cl.sol model.mps
I know that Gurobi does not necessarily approximate uniformly. In my previous post I was just using this as one possible approximation that I would have expected to be used. Using wolframalpha I found that a uniform grid on
N = ceil(20/(2*sqrt(2*ABS_ERROR)))
sections leads a desired maximal absolute error. So I would have expected Gurobi to select this or a similar approximation. At least it seems to be more or less straight forward to find an explicit (in this case even uniform) pwl approximation that fulfills an absolute error of at most 20. This was only for trying to understand what kind of appoximation gurobi might have used under the hood (but not really successfull in finding it in every case).
Thank you for looking into this. Best,
Stefan
----
main.cpp used to generate the above output and files:#include <iostream>
#include "gurobi_c++.h"
#include <vector>
#define max_x 20
#define target_x 4
int main() {
GRBEnv env = new GRBEnv(true);
env.start();
GRBModel model = GRBModel(env);
const double max_y = max_x * max_x / 2;
GRBVar x = model.addVar(0, max_x, 0, GRB_CONTINUOUS, "x");
GRBVar y = model.addVar(0, max_y, 0, GRB_CONTINUOUS, "y");
GRBVar z = model.addVar(0, GRB_INFINITY, 0, GRB_CONTINUOUS, "z");
GRBLinExpr obj = z;
model.setObjective(obj, GRB_MINIMIZE);
// model y = 1/2 * x^2 using a piecewise linear function
const double p[3]{0.5, 0, 0};
model.addGenConstrPoly(x, y, 3, p, "gen_const",
"FuncPieces=-1 FuncPieceError=20");
// Making sure that target_x is the optimal solution
model.addConstr(z, GRB_GREATER_EQUAL, target_x - x, "z_constr_1");
model.addConstr(z, GRB_GREATER_EQUAL, x - target_x, "z_constr_2");
model.optimize();
model.write("model.lp");
model.write("model.mps");
model.write("model.sol");
auto x_val = x.get(GRB_DoubleAttr_X);
auto y_val = y.get(GRB_DoubleAttr_X);
auto constraint_violation = model.get(GRB_DoubleAttr_ConstrVio);
auto real_val = x_val * x_val / 2;
auto gen_constr_approximation_error = std::abs(y_val - real_val);
std::cout << "x = " << x_val << std::endl;
std::cout << "y = " << y_val << " (exact value: " << real_val << ")" << std::endl;
std::cout << "approximation error at x is " << gen_constr_approximation_error << std::endl;
std::cout << "gurobi specifies a constraint violation of " << constraint_violation << std::endl;
}0 -
-
-
Hi Stefan,
Thank you for the model and your insights. I can reproduce the behavior on my end and will investigate.
I will come back once I have something relevant.
Best regards,
Jaromił0 -
-
-
Hi Stefan,
We identified the source of the issue. A wrong violation calculation is performed and reported. The actual violation is \(< 20\). This issue will be fixed an upcoming release. Unfortunately, there is currently no workaround for this.
Thank you for reporting this issue.
Best regards,
Jaromił0 -
Please sign in to leave a comment.
Comments
6 comments