I'm stumped. I have a problem formulated for NLOPT in R. The current problem solves for 180 variables with 28 equality constraints
The code is re-used from a simpler version of the problem, earlier in my script, with 36 variables and 20 equality constraints that solves instantly using NLOPT_LD_SLSQP
as the algorithm.
The larger version of the problem with 180 variables produces the following, immediately, when using NLOPT_LD_SLSQP
:
NLopt solver status: -4 ( NLOPT_ROUNDOFF_LIMITED: Roundoff errors led
to a breakdown of the optimization algorithm. In this case, the
returned minimum may
still be useful. (e.g. this error occurs in NEWUOA if one tries to
achieve a tolerance too close to machine precision.) )
This puzzled me, given that it worked on the smaller version of the problem. As well, it returns the starting values and doesn't actually complete any iterations. So I implemented, NLOPT_AUGLAG_LD_EQ
as the main algorithm, and NLOPT_LD_SLSQP
as the local algorithm. Now the problem fails to solve and produces this:
NLopt solver status: -1 ( NLOPT_FAILURE: Generic failure code. )
If I reduce the tolerance it just fails faster... I took the toy problem into Excel to see if I had failed to formulate properly or if it was somehow infeasible, but it solved immediately. I can give you this file if you want. I took the Excel solution values and populated the functions in R and sure enough my constraints, and objective function appear to be just fine.
I wonder if someone can help me with this problem. Here is some code that will produce the issue (tested and confirmed) for anyone in R:
library(pracma)
library(nloptr)
#My constraint function uses the following:
#RHS of the equality constraints
f.rhs <- c(590.0000,4781.0000,4414.0000,120.0000,224.0000,
849.0000,4693.0000,4374.0000,85.0000,697.0000,0.0000,
0.0000,0.0000,1092.0000,1434.0000,2251.0133,3482.9867,
2316.1813,1873.8187,1622.1450,1206.8550,1240.0000,1233.0000,
933.8532,733.1468,486.7907,395.2093,526.0000)
#matrix of constraints
#the first 10 rows are row total constraints
#the remaining 18 constraints are column total constraints
#this is a 28x180 matrix. It's sort of big to have it in this
#code window, but this code should produce the matrix for you
conmat <- c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.101385354477184,0,0,0,0,0,0,0,0,0,0,0,0,0.101385354477184,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.879602416076347,0,0,0,0,0,0,0,0,0,0,0,0,0,0.879602416076347,0,0,0,0,0,0,0,0,0,0,0,0,0,9.72187831641552,0,0,0,0,0,0,0,0,0,0,0,0,0,0,9.72187831641552,0,0,0,0,0,0,0,0,0,0,0,0,50.825662951981,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,50.825662951981,0,0,0,0,0,0,0,0,0,0,0,61.4161196898944,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,61.4161196898944,0,0,0,0,0,0,0,0,0,0,76.5722969856399,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,76.5722969856399,0,0,0,0,0,0,0,0,0,85.0874667314792,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,85.0874667314792,0,0,0,0,0,0,0,0,70.1228611430807,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,70.1228611430807,0,0,0,0,0,0,0,72.2969630445657,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,72.2969630445657,0,0,0,0,0,0,70.9315070452785,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,70.9315070452785,0,0,0,0,0,54.6520210670868,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,54.6520210670868,0,0,0,0,44.0086494626126,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,44.0086494626126,0,0,0,20.019587567467,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,20.019587567467,0,0,14.1724093345295,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,14.1724093345295,0,10.4922206705268,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,10.4922206705268,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,3.56379085643383,0,0,0,0,0,0,0,0,0,0,0,3.56379085643383,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,36.3843642437252,0,0,0,0,0,0,0,0,0,0,0,0,36.3843642437252,0,0,0,0,0,0,0,0,0,0,0,0,0,0,208.581934690648,0,0,0,0,0,0,0,0,0,0,0,0,0,208.581934690648,0,0,0,0,0,0,0,0,0,0,0,0,0,649.993541449925,0,0,0,0,0,0,0,0,0,0,0,0,0,0,649.993541449925,0,0,0,0,0,0,0,0,0,0,0,0,620.425879840303,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,620.425879840303,0,0,0,0,0,0,0,0,0,0,0,532.113517313307,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,532.113517313307,0,0,0,0,0,0,0,0,0,0,487.289086271457,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,487.289086271457,0,0,0,0,0,0,0,0,0,355.571461649492,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,355.571461649492,0,0,0,0,0,0,0,0,370.187737611463,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,370.187737611463,0,0,0,0,0,0,0,377.604110342457,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,377.604110342457,0,0,0,0,0,0,286.391885309974,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,286.391885309974,0,0,0,0,0,230.617639447808,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,230.617639447808,0,0,0,0,150.901768695456,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,150.901768695456,0,0,0,106.827457261503,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,106.827457261503,0,0,102.516716725934,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,102.516716725934,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,478.206023100627,0,0,0,0,0,0,0,0,0,0,478.206023100627,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,656.900445375174,0,0,0,0,0,0,0,0,0,0,0,656.900445375174,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,835.409246572897,0,0,0,0,0,0,0,0,0,0,0,0,835.409246572897,0,0,0,0,0,0,0,0,0,0,0,0,0,0,913.757767618374,0,0,0,0,0,0,0,0,0,0,0,0,0,913.757767618374,0,0,0,0,0,0,0,0,0,0,0,0,0,483.249080637494,0,0,0,0,0,0,0,0,0,0,0,0,0,0,483.249080637494,0,0,0,0,0,0,0,0,0,0,0,0,317.775206732195,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,317.775206732195,0,0,0,0,0,0,0,0,0,0,0,237.647569920133,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,237.647569920133,0,0,0,0,0,0,0,0,0,0,144.341667617821,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,144.341667617821,0,0,0,0,0,0,0,0,0,124.724396017039,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,124.724396017039,0,0,0,0,0,0,0,0,98.7489023828154,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,98.7489023828154,0,0,0,0,0,0,0,55.1727766234815,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,55.1727766234815,0,0,0,0,0,0,44.4279889177619,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,44.4279889177619,0,0,0,0,0,25.9420568421214,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,25.9420568421214,0,0,0,0,18.3650860591977,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,18.3650860591977,0,0,0,18.8065580202466,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,18.8065580202466,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.625528507812058,0,0,0,0,0,0,0,0,0,0.625528507812058,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1.14682340323878,0,0,0,0,0,0,0,0,0,0,1.14682340323878,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,8.19022941749738,0,0,0,0,0,0,0,0,0,0,0,8.19022941749738,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,25.5271393599384,0,0,0,0,0,0,0,0,0,0,0,0,25.5271393599384,0,0,0,0,0,0,0,0,0,0,0,0,0,0,30.1338658434231,0,0,0,0,0,0,0,0,0,0,0,0,0,30.1338658434231,0,0,0,0,0,0,0,0,0,0,0,0,0,33.325484734614,0,0,0,0,0,0,0,0,0,0,0,0,0,0,33.325484734614,0,0,0,0,0,0,0,0,0,0,0,0,34.6159582253122,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,34.6159582253122,0,0,0,0,0,0,0,0,0,0,0,25.9425263509155,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,25.9425263509155,0,0,0,0,0,0,0,0,0,0,28.6721543150432,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,28.6721543150432,0,0,0,0,0,0,0,0,0,24.6001930717219,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,24.6001930717219,0,0,0,0,0,0,0,0,19.4063381820094,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,19.4063381820094,0,0,0,0,0,0,0,15.6269927027328,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,15.6269927027328,0,0,0,0,0,0,7.86354283325046,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,7.86354283325046,0,0,0,0,0,5.56681537403581,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5.56681537403581,0,0,0,0,4.09245697782775,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,4.09245697782775,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.0478232804137659,0,0,0,0,0,0,0,0,0.0478232804137659,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.4097388469824,0,0,0,0,0,0,0,0,0,0.4097388469824,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1.18088960908956,0,0,0,0,0,0,0,0,0,0,1.18088960908956,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2.13050439017778,0,0,0,0,0,0,0,0,0,0,0,2.13050439017778,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2.88537630873373,0,0,0,0,0,0,0,0,0,0,0,0,2.88537630873373,0,0,0,0,0,0,0,0,0,0,0,0,0,0,4.47104432861129,0,0,0,0,0,0,0,0,0,0,0,0,0,4.47104432861129,0,0,0,0,0,0,0,0,0,0,0,0,0,6.7066648766379,0,0,0,0,0,0,0,0,0,0,0,0,0,0,6.7066648766379,0,0,0,0,0,0,0,0,0,0,0,0,7.21292740479352,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,7.21292740479352,0,0,0,0,0,0,0,0,0,0,0,13.2288618280058,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,13.2288618280058,0,0,0,0,0,0,0,0,0,0,17.7096180695658,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,17.7096180695658,0,0,0,0,0,0,0,0,0,32.1794159026695,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,32.1794159026695,0,0,0,0,0,0,0,0,25.9125391288607,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,25.9125391288607,0,0,0,0,0,0,0,39.9666473126856,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,39.9666473126856,0,0,0,0,0,0,28.293474255415,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,28.293474255415,0,0,0,0,0,71.410473393917,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,71.410473393917,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.182608374202251,0,0,0,0,0,0,0,0.182608374202251,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2.83538237187642,0,0,0,0,0,0,0,0,2.83538237187642,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,19.2609656330223,0,0,0,0,0,0,0,0,0,19.2609656330223,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,80.7233037904926,0,0,0,0,0,0,0,0,0,0,80.7233037904926,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,85.6477154102395,0,0,0,0,0,0,0,0,0,0,0,85.6477154102395,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,99.9589781355406,0,0,0,0,0,0,0,0,0,0,0,0,99.9589781355406,0,0,0,0,0,0,0,0,0,0,0,0,0,0,111.246861310206,0,0,0,0,0,0,0,0,0,0,0,0,0,111.246861310206,0,0,0,0,0,0,0,0,0,0,0,0,0,95.4680594823716,0,0,0,0,0,0,0,0,0,0,0,0,0,0,95.4680594823716,0,0,0,0,0,0,0,0,0,0,0,0,110.647972395614,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,110.647972395614,0,0,0,0,0,0,0,0,0,0,0,121.020177240488,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,121.020177240488,0,0,0,0,0,0,0,0,0,0,82.2865998682753,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,82.2865998682753,0,0,0,0,0,0,0,0,0,63.071831445587,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,63.071831445587,0,0,0,0,0,0,0,0,28.5766379506572,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,28.5766379506572,0,0,0,0,0,0,0,26.2025287609293,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,26.2025287609293,0,0,0,0,0,0,24.3533348857096,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,24.3533348857096,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,6.77378361061052,0,0,0,0,0,0,6.77378361061052,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,66.8658893692034,0,0,0,0,0,0,0,66.8658893692034,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,301.030835688858,0,0,0,0,0,0,0,0,301.030835688858,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,790.527770127371,0,0,0,0,0,0,0,0,0,790.527770127371,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,571.353574962507,0,0,0,0,0,0,0,0,0,0,571.353574962507,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,465.451417348464,0,0,0,0,0,0,0,0,0,0,0,465.451417348464,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,385.48085698754,0,0,0,0,0,0,0,0,0,0,0,0,385.48085698754,0,0,0,0,0,0,0,0,0,0,0,0,0,0,298.799798756715,0,0,0,0,0,0,0,0,0,0,0,0,0,298.799798756715,0,0,0,0,0,0,0,0,0,0,0,0,0,305.715022425598,0,0,0,0,0,0,0,0,0,0,0,0,0,0,305.715022425598,0,0,0,0,0,0,0,0,0,0,0,0,296.891932828816,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,296.891932828816,0,0,0,0,0,0,0,0,0,0,0,189.45276364351,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,189.45276364351,0,0,0,0,0,0,0,0,0,0,145.213592426377,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,145.213592426377,0,0,0,0,0,0,0,0,0,62.942694262135,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,62.942694262135,0,0,0,0,0,0,0,0,57.7134986817454,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,57.7134986817454,0,0,0,0,0,0,0,35.5644312539887,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,35.5644312539887,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,600.919741246423,0,0,0,0,0,600.919741246423,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,664.169162623053,0,0,0,0,0,0,664.169162623053,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,847.398221128337,0,0,0,0,0,0,0,847.398221128337,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,903.249824396504,0,0,0,0,0,0,0,0,903.249824396504,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,402.366866953328,0,0,0,0,0,0,0,0,0,402.366866953328,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,280.141092267195,0,0,0,0,0,0,0,0,0,0,280.141092267195,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,201.892549678482,0,0,0,0,0,0,0,0,0,0,0,201.892549678482,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,139.225117820461,0,0,0,0,0,0,0,0,0,0,0,0,139.225117820461,0,0,0,0,0,0,0,0,0,0,0,0,0,0,123.454735574923,0,0,0,0,0,0,0,0,0,0,0,0,0,123.454735574923,0,0,0,0,0,0,0,0,0,0,0,0,0,113.363999137925,0,0,0,0,0,0,0,0,0,0,0,0,0,0,113.363999137925,0,0,0,0,0,0,0,0,0,0,0,0,72.6983780508768,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,72.6983780508768,0,0,0,0,0,0,0,0,0,0,0,55.7225581581023,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,55.7225581581023,0,0,0,0,0,0,0,0,0,0,31.2411980593257,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,31.2411980593257,0,0,0,0,0,0,0,0,0,28.6457207488449,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,28.6457207488449,
0,0,0,0,0,0,0,0,38.4996441218375,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,38.4996441218375,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1.10058560667843,0,0,0,0,1.10058560667843,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,3.80905898924575,0,0,0,0,0,3.80905898924575,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,18.1504505832434,0,0,0,0,0,0,18.1504505832434,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,61.5743643896698,0,0,0,0,0,0,0,61.5743643896698,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,52.0421527000771,0,0,0,0,0,0,0,0,52.0421527000771,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,53.6771028451664,0,0,0,0,0,0,0,0,0,53.6771028451664,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,54.7413143789797,0,0,0,0,0,0,0,0,0,0,54.7413143789797,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,46.2299586575942,0,0,0,0,0,0,0,0,0,0,0,46.2299586575942,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,43.5728165091606,0,0,0,0,0,0,0,0,0,0,0,0,43.5728165091606,0,0,0,0,0,0,0,0,0,0,0,0,0,0,36.8659311062071,0,0,0,0,0,0,0,0,0,0,0,0,0,36.8659311062071,0,0,0,0,0,0,0,0,0,0,0,0,0,21.583566050924,0,0,0,0,0,0,0,0,0,0,0,0,0,0,21.583566050924,0,0,0,0,0,0,0,0,0,0,0,0,16.5435811193775,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,16.5435811193775,0,0,0,0,0,0,0,0,0,0,0,6.42202062218569,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,6.42202062218569,0,0,0,0,0,0,0,0,0,0,5.88848766417714,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5.88848766417714,0,0,0,0,0,0,0,0,0,3.34484908130525,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,3.34484908130525,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.478730062097792,0,0,0,0.478730062097792,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.599532361620929,0,0,0,0,0.599532361620929,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2.08861096083086,0,0,0,0,0,2.08861096083086,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,4.67685892604545,0,0,0,0,0,0,4.67685892604545,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,6.66062843004317,0,0,0,0,0,0,0,6.66062843004317,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,10.3325985324682,0,0,0,0,0,0,0,0,10.3325985324682,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,17.4366302473174,0,0,0,0,0,0,0,0,0,17.4366302473174,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,23.9406624891688,0,0,0,0,0,0,0,0,0,0,23.9406624891688,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,47.4993402782371,0,0,0,0,0,0,0,0,0,0,0,47.4993402782371,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,75.2636287745429,0,0,0,0,0,0,0,0,0,0,0,0,75.2636287745429,0,0,0,0,0,0,0,0,0,0,0,0,0,0,120.029503824241,0,0,0,0,0,0,0,0,0,0,0,0,0,120.029503824241,0,0,0,0,0,0,0,0,0,0,0,0,0,92.0013786669861,0,0,0,0,0,0,0,0,0,0,0,0,0,0,92.0013786669861,0,0,0,0,0,0,0,0,0,0,0,0,112.914580678886,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,112.914580678886,0,0,0,0,0,0,0,0,0,0,0,103.53378703525,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,103.53378703525,0,0,0,0,0,0,0,0,0,0,216.919314868552,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,216.919314868552)
#Create the matrix from my list of values
conmat <- matrix(conmat,nrow=28,ncol=180)
#Create the constraint function so that it produces 0s
#for the equality constraints
eqn <- function(x){
z=c()
for (i in 1:length(f.rhs)){
z[i]=x%*%conmat[i,]-f.rhs[i]
}
return(z)
}
#Function for the Jacobian of the constraint function
eqn_grad <- function(x){
jacobian(eqn,x)
}
#Create my objective function: sum of squared error
#Data for the function is in f.obj
f.obj<- c(0,0,0,0.101385354477184,0.879602416076347,9.72187831641552,50.825662951981,61.4161196898944,76.5722969856399,85.0874667314792,70.1228611430807,72.2969630445657,70.9315070452785,54.6520210670868,44.0086494626126,20.019587567467,14.1724093345295,10.4922206705268,0,0,0,3.56379085643383,36.3843642437252,208.581934690648,649.993541449925,620.425879840303,532.113517313307,487.289086271457,355.571461649492,370.187737611463,377.604110342457,286.391885309974,230.617639447808,150.901768695456,106.827457261503,102.516716725934,0,0,0,478.206023100627,656.900445375174,835.409246572897,913.757767618374,483.249080637494,317.775206732195,237.647569920133,144.341667617821,124.724396017039,98.7489023828154,55.1727766234815,44.4279889177619,25.9420568421214,18.3650860591977,18.8065580202466,0,0,0,0.625528507812058,1.14682340323878,8.19022941749738,25.5271393599384,30.1338658434231,33.325484734614,34.6159582253122,25.9425263509155,28.6721543150432,24.6001930717219,19.4063381820094,15.6269927027328,7.86354283325046,5.56681537403581,4.09245697782775,0,0,0,0.0478232804137659,0.4097388469824,1.18088960908956,2.13050439017778,2.88537630873373,4.47104432861129,6.7066648766379,7.21292740479352,13.2288618280058,17.7096180695658,32.1794159026695,25.9125391288607,39.9666473126856,28.293474255415,71.410473393917,0,0,0,0.182608374202251,2.83538237187642,19.2609656330223,80.7233037904926,85.6477154102395,99.9589781355406,111.246861310206,95.4680594823716,110.647972395614,121.020177240488,82.2865998682753,63.071831445587,28.5766379506572,26.2025287609293,24.3533348857096,0,0,0,6.77378361061052,66.8658893692034,301.030835688858,790.527770127371,571.353574962507,465.451417348464,385.48085698754,298.799798756715,305.715022425598,296.891932828816,189.45276364351,145.213592426377,62.942694262135,57.7134986817454,35.5644312539887,0,0,0,600.919741246423,664.169162623053,847.398221128337,903.249824396504,402.366866953328,280.141092267195,201.892549678482,139.225117820461,123.454735574923,113.363999137925,72.6983780508768,55.7225581581023,31.2411980593257,28.6457207488449,38.4996441218375,0,0,0,1.10058560667843,3.80905898924575,18.1504505832434,61.5743643896698,52.0421527000771,53.6771028451664,54.7413143789797,46.2299586575942,43.5728165091606,36.8659311062071,21.583566050924,16.5435811193775,6.42202062218569,5.88848766417714,3.34484908130525,0,0,0,0.478730062097792,0.599532361620929,2.08861096083086,4.67685892604545,6.66062843004317,10.3325985324682,17.4366302473174,23.9406624891688,47.4993402782371,75.2636287745429,120.029503824241,92.0013786669861,112.914580678886,103.53378703525,216.919314868552)
#Objective Function is SSE
fn <- function(x){
sum((f.obj*x-f.obj)^2)
}
#Objective Gradient
fn_grad <- function(x){
grad(fn,x)
}
#Optimization options:
#Starting Values
x0 <- c(matrix(1,1,ncol=length(f.obj)))
#Lower Bound
lb_x <- c(matrix(0,1,ncol=length(f.obj)))
#Upper Bound
ub_x <- c(matrix(2,1,ncol=length(f.obj)))
Now here is the code to solve with SLSQP:
opts_list <- list("algorithm"="NLOPT_LD_SLSQP", "xtol_rel"=1.0e-8,"maxeval"=1000)
SOLUTION <- nloptr(x0,
eval_f=fn,
eval_grad_f=fn_grad,
lb=lb_x,
ub=ub_x,
eval_g_eq = eqn,
eval_jac_g_eq = eqn_grad,
opts=opts_list)
SOLUTION
And here is the code to solve with AUGLAG:
local_opts_list <- list("algorithm" = "NLOPT_LD_SLSQP","xtol_rel"=1.0e-8)
opts_list <- list("algorithm"="NLOPT_LD_AUGLAG_EQ", "xtol_rel"=1.0e-8,"maxeval"=1000, "local_opts"=local_opts_list)
SOLUTION <- nloptr(x0,
eval_f=fn,
eval_grad_f=fn_grad,
lb=lb_x,
ub=ub_x,
eval_g_eq = eqn,
eval_jac_g_eq = eqn_grad,
opts=opts_list)
SOLUTION
Last, but not least, here is a feasible solution produced by Excel:
XL_Sol <- c(1,1,1,0.999840758,0.998619343,0.984736297,0.917305327,0.895766204,0.902054225,0.930449654,0.908095562,0.90555472,0.907109309,0.92775884,0.941019433,0.973119106,0.980882453,1.00818501,1,1,1,0.999277369,0.99265595,0.957137798,1.068940265,0.971466259,0.916751383,1.031659407,1.088943875,1.106986398,1.09298138,1.169120549,1.136193658,1.242584994,1.260102262,1.218743178,1,1,1,1.036812057,1.051986714,1.066969218,0.764543983,0.983061755,1.117743946,1.086522598,1.048114102,1.042021839,1.032920906,1.017781177,1.013537363,1.007869529,1.005459785,1.04562595,1,1,1,0.986367866,0.975008481,0.821522583,0.442429864,0.339768188,0.28380776,0.279905133,0.44151851,0.382930646,0.470415572,0.581960382,0.663070622,0.830426413,0.879918273,0.920434632,1,1,1,0.99976331,0.997972512,0.994156123,0.989331379,0.985345491,0.979153535,0.96871128,0.966148601,0.937972897,0.916916787,0.848764035,0.877734945,0.811485291,0.866305968,0.895221739,1,1,1,0.999728838,0.995792558,0.971414471,0.875797499,0.862210736,0.88101002,0.867466136,0.883358551,0.865529128,0.894845599,0.898520124,0.920974341,0.964087977,0.966910775,1.02107896,1,1,1,1.000481883,1.004812176,0.752446473,1.516261232,1.234618416,1.169698414,1.138623218,1.098571557,1.102118823,1.098211297,1.060360506,1.043658294,1.01881129,1.016899233,1.086075848,1,1,1,0.970757235,0.948975599,1.042336122,0.836515328,0.953263454,1.082324289,1.058181611,1.036019128,1.032395622,1.02937951,1.018031769,1.012837354,1.007153339,1.006385936,1.090482373,1,1,1,0.979339786,0.92850077,0.659360317,0,0.01800471,0.038870008,0,0.145186371,0.194368327,0.318022546,0.600269991,0.693283542,0.880910409,0.890768254,0.94507392,1,1,1,0.998868123,0.998583119,0.995062994,0.988669793,0.983388904,0.978533683,0.963729915,0.94954262,0.90018972,0.84196737,0.764585537,0.80459333,0.760693136,0.77966356,0.896573068)
#eqn(XL_Sol) should give you a vector of numbers that are pretty close to zero.
My question is: why does this model produce the errors I've reported here for each of the SLSQP, and AUGLAG algorithms implemented in the sample code?
Would love some input here. Please let me know if you need any additional information!
An interesting thing about your system is that conmat
is rank deficient:
qr(conmat)$rank # 24 < 28 (number of rows in the system)
But the system conmat%*%x = f.rhs
is actually consistent and you can find one of its solutions using the answers from this post.
That means that your constraint system has 4 redundant relations. Sometimes issues like this can lead to numerical instability, so let's see what happens if we remove the redundant equations. To do this, we can adapt the answer in this post:
# let's relabel
conmat_raw = conmat
f.rhs_raw = f.rhs
# append rhs as a column to conmat
augmented_matrix <- cbind(conmat, f.rhs)
# qr transform of the transpose
q <- qr(t(augmented_matrix))
# extract a set of linearly dependent rows for the original augmented matrix
reduced <- t(t(augmented_matrix)[,q$pivot[seq(q$rank)]])
# decompose augmented matrix
conmat = reduced[,1:180]
f.rhs = reduced[,181]
Now, we can run the NLOPT_LD_SLSQP
solver. After 1350 iterations, I get a solution. Let's see how it compares to the XL_Sol
.
> sol = SOLUTION$solution
> fn(sol)
[1] 58404.02
> fn(XL_Sol)
[1] 317837
So we've improved against the objective. Let's check the solutions using the original constraints:
> sum((conmat_raw%*%sol - f.rhs_raw)^2)
[1] 6.071796e-08
> sum((conmat_raw%*%XL_Sol - f.rhs_raw)^2)
[1] 2.771957e-08
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With