The return value doesnt satisfy my inequality constraints (I set the precision
to 1e-8, and the solution evaluates the inequality constraint to +5e-7)
In this particular case, it returns after reaching maxeval (it seems to
oscillate infinitely near the constraint boundary), but I think I saw that
happen with "xtol reached" as well. The attached code illustrates that.
thanks
Christophe
-----Original Message-----
From: [email protected]
[mailto:[email protected]] On Behalf Of Steven G. Johnson
Sent: 21 November 2011 21:59
To: nlopt-discuss
Subject: Re: [NLopt-discuss] MMA method moving out of feasible region
On Nov 21, 2011, at 9:06 AM, chris wrote:
> I'm using NLOpt to solve a (strictly convex) QP problem with box,
> linear and
> quadratic inequality constraints. I tried the MMA and SLSQP methods.
> In most
> cases, they both converge to a very good solution but SLSQP always
> gives an
> answer that satisfies the inequality constraints up to required
> precision,
> whereas MMA occasionally starts erring into the infeasible region
> (despite
> starting with a feasible point). Did I misunderstand the strict
> requirement of
> feasibility of NLOpt ? (I can send an example .c file if that helps).
In general, NLopt only guarantees strict feasibility of all
evaluations for box constraints (specified explicitly as lower/upper
bounds); with other constraints, you should get something that
converges towards a feasible point (to the attainable accuracy), but
may stray outside the feasible region at intermediate steps.
MMA, in particular, has both "outer iterations" and "inner iterations"
within each outer iteration. The outer iterations are guaranteed to
be feasible (from a feasible starting point), but the inner iterations
may evaluate the function at points outside the feasible region. It
should only return the result of an outer iteration, however, i.e. a
feasible point.
When you say "MMA occasionally starts erring into the infeasible
region" are you talking about the return value of the optimization, or
just about intermediate steps prior to return? In the latter case
there is no guarantee of feasibility.
Steven
_______________________________________________
NLopt-discuss mailing list
[email protected]
http://ab-initio.mit.edu/cgi-bin/mailman/listinfo/nlopt-discuss
#include <vector>
#include <nlopt.hpp>
#include <boost/numeric/ublas/symmetric.hpp>
using std::vector;
namespace ublas = boost::numeric::ublas;
//data
double c_initialGuess [10] =
{0.00000000000000000,3.3751120535445790e-005,7.8331354472156156e-006,3.7022683684448233e-005,0.00000000000000000,0.00000000000000000,0.00000000000000000,0.00000000000000000,0.00000000000000000,0.99992139306033290};
double c_Returns [10] =
{0.57534183730894939,0.75361853693561376,0.81320707241841605,0.68584281649105239,1.0867187892134469,0.73875081197658521,1.4198848137587050,0.94029934796451631,1.6158755405308778,0.0017001496300403673};
double c_PreviousWeights [10] =
{7.2310069527092565e-010,0.081617502360159688,9.7299353070326459e-010,9.6761777672132411e-010,0.093434565210364390,0.035697742262116167,0.11497477392474496,0.0027518721039319952,0.066710659656591662,0.60481288181837911};
double c_TradingCosts [10] =
{0.00000000000000000,0.00000000000000000,0.00000000000000000,0.00000000000000000,0.00000000000000000,0.00000000000000000,0.00000000000000000,0.00000000000000000,0.00000000000000000,0.00000000000000000};
double c_LinearConstraints [10] =
{0.24413374983701874,0.23751704919147126,0.27536692441768351,0.22589935478529985,0.25980419238619457,0.23935081313915693,0.19673766819835123,0.24464901064749398,0.32654123739732954,0.00000000000000000};
double c_MinWeights [10] =
{0.00000000000000000,0.00000000000000000,0.00000000000000000,0.00000000000000000,0.00000000000000000,0.00000000000000000,0.00000000000000000,0.00000000000000000,0.00000000000000000,0.52354384150084243};
double c_MaxWeights [10] =
{0.094885101320799775,0.10304880883975288,0.079046830901594312,0.098721896375702167,0.092032500019005037,0.093054552061202084,0.11434947803979781,0.093898207342585319,0.066683627024280934,1.0000000000000000};
double c_maxVol = 0.13;
double c_linearMax = 0.11911403962478939;
double c_CovarMat [55] =
{0.171229891970579,0.156300524872849,0.160990718302128,0.138458782868199,0.115986075282479,0.216821131395112,0.026003182237491,0.0287274799318514,0.0434614126903358,0.145810329295916,0.06499315086136,0.0566178116838743,0.110069964275089,0.0961393900426599,0.190512091047068,0.0261232948421907,0.032368123956648,0.0455856170670122,0.0396987740240037,0.116026900129326,0.1636173408517,0.0929812467599921,0.0838528635981951,0.105241713004538,0.0706972347447946,0.0886511861868356,0.0435217368285407,0.103687318436377,0.0928591568313496,0.0814136377069847,0.130814061606889,0.0717113696849801,0.155245392269717,0.0922462030623368,0.100178247097967,0.169617682356235,0.118778573528206,0.0987988496099389,0.165512982670255,0.095543125901775,0.154449330899773,0.109798081628768,0.110266855359186,0.154274703455513,0.297844906933627,-6.38916323416339E-06,-7.40291721073104E-06,-7.21929002945255E-06,-6.70545567630249E-06,-3.5693585953525E-06,2.18812778817298E-06,-5.56016476978313E-06,-4.60194495732275E-06,-3.82704986022117E-06,3.93985877793975E-09};
class GPSOptimization
{
public:
GPSOptimization(int n,
double* covarMat,
double* returns,
double* previousWeights,
double* tradingCosts,
double* linearConstraint,
double* minWeights,
double* maxWeights,
double linearMax,
double
maxVol):myReturns(n),myPreviousWeights(n),myLinearConstraints(n),myTradingCosts(n),myMinWeights(n),myMaxWeights(n),myCovarMat(n)
{
size_t k=0;
for(size_t i=0; i<n; i++)
for(size_t j=0; j<=i; j++)
{
myCovarMat(i,j)=covarMat[k];
k++;
}
for(size_t i=0; i<n; i++)
{
myReturns[i] = returns[i];
myPreviousWeights[i]=previousWeights[i];
myTradingCosts[i]=tradingCosts[i];
myLinearConstraints[i]=linearConstraint[i];
myMinWeights[i]=minWeights[i];
myMaxWeights[i]=maxWeights[i];
}
myLinearMax = linearMax;
myMaxVol = maxVol;
}
const vector<double>& getReturns() const {return myReturns;}
const vector<double>& getPreviousWeights() const {return
myPreviousWeights;}
const vector<double>& getTradingCosts() const {return myTradingCosts;}
const vector<double>& getLinearConstraints() const {return
myLinearConstraints;}
const vector<double>& getMinWeights() const {return myMinWeights;}
const vector<double>& getMaxWeights() const {return myMaxWeights;}
double getLinearMax() const {return myLinearMax;}
double getMaxVol() const {return myMaxVol;}
const ublas::symmetric_matrix<double>& getCovarianceMatrix() const
{return myCovarMat;}
private:
vector<double>
myReturns,myPreviousWeights,myTradingCosts,myLinearConstraints,myMinWeights,myMaxWeights;
double myMaxVol,myLinearMax;
ublas::symmetric_matrix<double> myCovarMat;
};
//follow NLOpt type
double NLOptObjectiveFunction(size_t n, const double* x, double* grad, void*
f_data)
{
GPSOptimization* me = (GPSOptimization*)(f_data);
const std::vector<double>& returns = me->getReturns();
const std::vector<double>& tradingCosts = me->getTradingCosts();
const std::vector<double>& previousWeights = me->getPreviousWeights();
double res=0, sumx=0;
for(size_t i=0; i<n; ++i)
{
sumx += x[i];
res += returns[i]*x[i];
res -=
(x[i]-previousWeights[i])*(x[i]-previousWeights[i])*tradingCosts[i];
if(grad)
{
grad[i]=returns[i]-2*(x[i]-previousWeights[i])*tradingCosts[i];
}
}
double x_n = 1-sumx;
res += returns[n]*x_n;
res -=
(x_n-previousWeights[n])*(x_n-previousWeights[n])*tradingCosts[n];
if(grad)
for(size_t i=0; i<n; ++i)
grad[i] +=
-returns[n]+2*(x_n-previousWeights[n])*tradingCosts[n];
return res;
}
double NLOptVolatilityInequalityConstraints(size_t n, const double* x, double*
grad, void* f_data)
{
GPSOptimization* me = (GPSOptimization*)(f_data);
const ublas::matrix<double>& covarianceMat = me->getCovarianceMatrix();
double maxVol = me->getMaxVol();
ublas::vector<double> X(n+1);
X[n]=1;
for (size_t i=0; i<n; i++)
{
X[i]=x[i];
X[n]-=x[i];
}
ublas::vector<double> Y = ublas::prod(covarianceMat,X);
if(grad)
for(size_t i=0; i<n; i++)
grad[i]=2*(Y[i]-Y[n]);
double res=0;
for(size_t i=0; i<n; i++)
res += X[i]*Y[i];
res += X[n]*Y[n];
res-=maxVol*maxVol;
return res;
}
double NLOptLinearInequalityConstraints(size_t n, const double* x, double*
grad, void* f_data)
{
GPSOptimization* me = (GPSOptimization*)(f_data);
const std::vector<double>& linearConstraint =
me->getLinearConstraints();
double linearMax = me->getLinearMax();
double res =linearConstraint[n];
for (size_t i=0; i<n; i++)
{
res += x[i]*(linearConstraint[i]-linearConstraint[n]);
if(grad)
grad[i]=(linearConstraint[i]-linearConstraint[n]);
}
res-=linearMax;
return res;
}
double NLOptLinearEqualityConstraintsUp(size_t n, const double* x, double*
grad, void* f_data)
{
GPSOptimization* me = (GPSOptimization*)(f_data);
double floor = me->getMinWeights().back();
double res =0;
for (size_t i=0; i<n; i++)
{
res += x[i];
if(grad)
grad[i]=1.;
}
res+=floor-1;
return res;
}
double NLOptLinearEqualityConstraintsDown(size_t n, const double* x, double*
grad, void* f_data)
{
GPSOptimization* me = (GPSOptimization*)(f_data);
double cap = me->getMaxWeights().back();
double res =0;
for (size_t i=0; i<n; i++)
{
res -= x[i];
if(grad)
grad[i]=-1.;
}
res += 1-cap;
return res;
}
bool computeNLOpt(vector<double>& weights, const vector<double>& initialGuess,
const GPSOptimization* gps, double epsilon)
{
size_t n = weights.size()-1;//parametrize x_n as 1-sum(x_i) to remove
equality constraint
nlopt::opt optim(nlopt::LD_MMA,n);
if(NLOptVolatilityInequalityConstraints(n,&initialGuess[0],NULL,const_cast<GPSOptimization*>(gps))>epsilon)
throw "Constraint violation in initial guess";
if(NLOptLinearInequalityConstraints(n,&initialGuess[0],NULL,const_cast<GPSOptimization*>(gps))>epsilon)
"Constraint violation in initial guess";
if(NLOptLinearEqualityConstraintsUp(n,&initialGuess[0],NULL,const_cast<GPSOptimization*>(gps))>epsilon)
"Constraint violation in initial guess";
if(NLOptLinearEqualityConstraintsDown(n,&initialGuess[0],NULL,const_cast<GPSOptimization*>(gps))>epsilon)
"Constraint violation in initial guess";
optim.set_max_objective(&NLOptObjectiveFunction,const_cast<GPSOptimization*>(gps));
optim.add_inequality_constraint(&NLOptVolatilityInequalityConstraints,const_cast<GPSOptimization*>(gps),epsilon);
optim.add_inequality_constraint(&NLOptLinearInequalityConstraints,const_cast<GPSOptimization*>(gps),epsilon);
optim.add_inequality_constraint(&NLOptLinearEqualityConstraintsUp,const_cast<GPSOptimization*>(gps),epsilon);
optim.add_inequality_constraint(&NLOptLinearEqualityConstraintsDown,const_cast<GPSOptimization*>(gps),epsilon);
optim.set_lower_bounds(std::vector<double>(gps->getMinWeights().begin(),gps->getMinWeights().begin()+n));
optim.set_upper_bounds(std::vector<double>(gps->getMaxWeights().begin(),gps->getMaxWeights().begin()+n));
optim.set_xtol_abs(epsilon);
optim.set_maxeval(1000);
//optim.set_maxtime(10.0);
weights =
std::vector<double>(initialGuess.begin(),initialGuess.begin()+n);
double optimizedValue=0;
nlopt::result res = optim.optimize(weights,optimizedValue);
//It will sometimes return a result even when the constraints cannot be
matched. We test for that.
if(NLOptVolatilityInequalityConstraints(n,&weights[0],NULL,const_cast<GPSOptimization*>(gps))>epsilon)
throw "Constraint violation in result";
if(NLOptLinearInequalityConstraints(n,&weights[0],NULL,const_cast<GPSOptimization*>(gps))>epsilon)
"Constraint violation in result";
if(NLOptLinearEqualityConstraintsUp(n,&weights[0],NULL,const_cast<GPSOptimization*>(gps))>epsilon)
"Constraint violation in result";
if(NLOptLinearEqualityConstraintsDown(n,&weights[0],NULL,const_cast<GPSOptimization*>(gps))>epsilon)
"Constraint violation in result";
weights.push_back(1.);
for(size_t i=0; i<n; i++)
weights.back()-=weights[i];
return (res>0 && res<5); //success without maxeval or maxtime being
reached
}
int main()
{
GPSOptimization problem(10,c_CovarMat,
c_Returns,
c_PreviousWeights,
c_TradingCosts,
c_LinearConstraints,
c_MinWeights,
c_MaxWeights,
c_linearMax,
c_maxVol);
vector<double> initialGuess(c_initialGuess,c_initialGuess+10);
vector<double> weights = initialGuess;
try
{
computeNLOpt(weights,initialGuess,&problem,1e-8);
}
catch (char* e)
{
std::cout << e;
return 0;
}
for(size_t i=0; i<10; i++)
std::cout << weights[i] << "\n";
return 0;
}
_______________________________________________
NLopt-discuss mailing list
[email protected]
http://ab-initio.mit.edu/cgi-bin/mailman/listinfo/nlopt-discuss
