Program Listing for File qp_param.cpp
↰ Return to documentation for file (src/mathopt/mp_param/qp_param.cpp)
/* #############################################
* This file is part of
* ZERO
*
* Copyright (c) 2020
* Released under the Creative Commons
* CC BY-NC-SA 4.0 License
*
* Find out more at
* https://github.com/ds4dm/ZERO
* #############################################*/
#include "mathopt/mp_param/qp_param.h"
std::ostream &MathOpt::operator<<(std::ostream &os, const MathOpt::QP_Param &Q) {
os << "Quadratic program with linear inequality constraints: " << '\n';
os << Q.getNumVars() << " decision variables parametrized by " << Q.getNumParams() << " variables"
<< '\n';
os << Q.getb().n_rows << " linear inequalities" << '\n' << '\n';
return os;
}
bool MathOpt::QP_Param::operator==(const QP_Param &Q2) const {
if (!Utils::isZero(this->Q - Q2.getQ()))
return false;
if (!Utils::isZero(this->C - Q2.getC()))
return false;
if (!Utils::isZero(this->getA(true) - Q2.getA(true)))
return false;
if (!Utils::isZero(this->getB(true) - Q2.getB(true)))
return false;
if (!Utils::isZero(this->c - Q2.getc()))
return false;
if (!Utils::isZero(this->d - Q2.getd()))
return false;
for (unsigned int i = 0; i < this->Bounds.size(); ++i)
if (this->Bounds.at(i) != Q2.Bounds.at(i))
return false;
if (!Utils::isZero(this->getb(true) - Q2.getb(true)))
return false;
return true;
}
void MathOpt::QP_Param::makeyQy() {
if (this->MadeyQy)
return;
GRBVar y[this->numVars];
for (unsigned int i = 0; i < numVars; i++)
y[i] = this->Model.addVar(
Bounds.at(i).first, Bounds.at(i).second, 0, GRB_CONTINUOUS, "y_" + std::to_string(i));
GRBQuadExpr yQy{0};
for (auto val = Q.begin(); val != Q.end(); ++val) {
unsigned int i, j;
double value = (*val);
i = val.row();
j = val.col();
yQy += 0.5 * y[i] * value * y[j];
}
Model.setObjective(yQy, GRB_MINIMIZE);
Model.update();
this->MadeyQy = true;
}
std::unique_ptr<GRBModel> MathOpt::QP_Param::solveFixed(arma::vec x, bool solve) {
this->makeyQy();
if (x.size() != this->numParams)
throw ZEROException(ZEROErrorCode::Assertion,
"Mismatch in x size: " + std::to_string(x.size()) +
" != " + std::to_string(numParams));
std::unique_ptr<GRBModel> model(new GRBModel(this->Model));
try {
GRBQuadExpr yQy = model->getObjective()+ arma::as_scalar(this->d.t() * x);
arma::vec Cx, Ax;
Cx = this->C * x;
Ax = this->A * x;
GRBVar y[this->numVars];
for (unsigned int i = 0; i < this->numVars; i++) {
y[i] = model->getVarByName("y_" + std::to_string(i));
yQy += (Cx[i] + c[i]) * y[i];
}
model->setObjective(yQy , GRB_MINIMIZE);
Utils::addSparseConstraints(B, b - Ax, y, "Constr_", model.get(), GRB_LESS_EQUAL, nullptr);
model->update();
model->set(GRB_IntParam_OutputFlag, 0);
model->set(GRB_IntParam_NonConvex, 2);
if (solve)
model->optimize();
} catch (GRBException &e) {
throw ZEROException(e);
}
return model;
}
unsigned int MathOpt::QP_Param::KKT(arma::sp_mat &M, arma::sp_mat &N, arma::vec &q) const {
this->forceDataCheck();
M = arma::join_cols( // In armadillo join_cols(A, B) is same as [A;B] in
// Matlab
// join_rows(A, B) is same as [A B] in Matlab
arma::join_rows(this->Q, this->getB(true).t()),
arma::join_rows(-this->getB(true),
arma::zeros<arma::sp_mat>(this->numConstr, this->numConstr)));
ZEROAssert(M.n_cols == (numVars + numConstr + this->B_bounds.n_rows));
N = arma::join_cols(this->C, -this->getA(true));
ZEROAssert(N.n_cols == numParams);
q = arma::join_cols(this->c, this->getb(true));
ZEROAssert(q.size() == (this->c.size() + this->b.size() + this->b_bounds.size()));
// q.print();
return M.n_rows;
}
MathOpt::QP_Param &MathOpt::QP_Param::set(const arma::sp_mat &Q_in,
const arma::sp_mat &C_in,
const arma::sp_mat &A_in,
const arma::sp_mat &B_in,
const arma::vec &c_in,
const arma::vec &b_in,
const arma::vec &d_in = {}) {
this->MadeyQy = false;
MP_Param::set(Q_in, C_in, A_in, B_in, c_in, b_in, d_in);
return *this;
}
MathOpt::QP_Param &MathOpt::QP_Param::set(arma::sp_mat &&Q_in,
arma::sp_mat &&C_in,
arma::sp_mat &&A_in,
arma::sp_mat &&B_in,
arma::vec &&c_in,
arma::vec &&b_in,
arma::vec &&d_in = {}) {
this->MadeyQy = false;
MP_Param::set(Q_in, C_in, A_in, B_in, c_in, b_in, d_in);
return *this;
}
MathOpt::QP_Param &MathOpt::QP_Param::set(QP_Objective &&obj, QP_Constraints &&cons) {
return this->set(std::move(obj.Q),
std::move(obj.C),
std::move(cons.A),
std::move(cons.B),
std::move(obj.c),
std::move(cons.b),
std::move(obj.d));
}
MathOpt::QP_Param &MathOpt::QP_Param::set(const QP_Objective &obj, const QP_Constraints &cons) {
return this->set(obj.Q, obj.C, cons.A, cons.B, obj.c, cons.b, obj.d);
}
void MathOpt::QP_Param::save(const std::string &filename, bool append) const {
Utils::appendSave(std::string("QP_Param"), filename, append);
Utils::appendSave(this->Q, filename, std::string("QP_Param::Q"), false);
Utils::appendSave(this->getA(true), filename, std::string("QP_Param::A"), false);
Utils::appendSave(this->getB(true), filename, std::string("QP_Param::B"), false);
Utils::appendSave(this->C, filename, std::string("QP_Param::C"), false);
Utils::appendSave(this->getb(true), filename, std::string("QP_Param::b"), false);
Utils::appendSave(this->c, filename, std::string("QP_Param::c"), false);
Utils::appendSave(this->d, filename, std::string("QP_Param::d"), false);
arma::sp_mat BO(this->numVars, 2);
for (unsigned int i = 0; i < this->numVars; ++i) {
BO.at(i, 0) = this->Bounds.at(i).first;
BO.at(i, 1) = this->Bounds.at(i).second;
}
Utils::appendSave(BO, filename, std::string("QP_Param::Bounds"), false);
LOG_S(1) << "Saved QP_Param to file " << filename;
}
long int MathOpt::QP_Param::load(const std::string &filename, long int pos) {
arma::sp_mat Q_in, A_in, B_in, C_in, BO;
arma::vec c_in, b_in, d_in;
std::string headercheck;
pos = Utils::appendRead(headercheck, filename, pos);
if (headercheck != "QP_Param")
throw ZEROException(ZEROErrorCode::IOError, "Invalid header");
pos = Utils::appendRead(Q_in, filename, pos, std::string("QP_Param::Q"));
pos = Utils::appendRead(A_in, filename, pos, std::string("QP_Param::A"));
pos = Utils::appendRead(B_in, filename, pos, std::string("QP_Param::B"));
pos = Utils::appendRead(C_in, filename, pos, std::string("QP_Param::C"));
pos = Utils::appendRead(b_in, filename, pos, std::string("QP_Param::b"));
pos = Utils::appendRead(c_in, filename, pos, std::string("QP_Param::c"));
pos = Utils::appendRead(d_in, filename, pos, std::string("QP_Param::d"));
pos = Utils::appendRead(BO, filename, pos, std::string("QP_Param::Bounds"));
if (BO.n_rows > 0) {
ZEROAssert(BO.n_cols == 2);
for (unsigned int i = 0; i < B_in.n_cols; ++i)
this->Bounds.push_back(
{abs(BO.at(i, 0)) < 1e20 ? BO.at(i, 0) : Utils::getSign(BO.at(i, 0)) * GRB_INFINITY,
abs(BO.at(i, 1)) < 1e20 ? BO.at(i, 1) : Utils::getSign(BO.at(i, 1)) * GRB_INFINITY});
int diff = B_in.n_cols - BO.n_rows;
for (unsigned int i = 0; i < diff; ++i)
this->Bounds.push_back({0, GRB_INFINITY});
}
LOG_S(1) << "Loaded QP_Param to file " << filename;
this->set(Q_in, C_in, A_in, B_in, c_in, b_in, d_in);
return pos;
}
MathOpt::QP_Param::QP_Param(const arma::sp_mat &Q_in,
const arma::sp_mat &C_in,
const arma::sp_mat &A_in,
const arma::sp_mat &B_in,
const arma::vec &c_in,
const arma::vec &b_in,
const arma::vec &d_in,
GRBEnv *env)
: MP_Param(env), MadeyQy{false}, Model{(*env)} {
this->MadeyQy = false;
this->set(Q_in, C_in, A_in, B_in, c_in, b_in, d_in);
this->size();
this->forceDataCheck();
}