NLsolve.jl icon indicating copy to clipboard operation
NLsolve.jl copied to clipboard

Function straying out of bounds in NLsolve

Open davidzentlermunro opened this issue 10 years ago • 6 comments

I'm getting a domain error when running the NLsolve code given below: "DomainError in ^ at math.jl:252" I think this is because mcpsolve is straying out of the zero one bounds I have given it, and so tries to evaluate the root of a negative number. Is this something that could/should be fixed? I may well be wrong in my diagnosis though - happy to be corrected.

David

using Distributions
using Devectorize
using Distances
using StatsBase
using NumericExtensions
using NLsolve

beta = 0.95;                                                               
lambda0 = .90;                                                             
lambda1 = 0.05;                                                           
mu = 2;                                                                     
rho = 5.56;                                                                 
xmin= 0.73;                                                                  
xmax = xmin+1;                                                             
z0 = 0.0;                                                                   
nu = 0.64;                                                                  
sigma = 0.023;                                                             
alpha = 2;                                                                  
TFP = 1;                                                                    
eta = 0.3;                                                                  
delta = 0.01;                                                                                                                          
amaximum=500;                                                              
mwruns=15;
gamma=0.5                                                                   
kappa = 1                                                                   
psi=0.5
prod=linspace(xmin,xmax,ns);                                                    
l1=0.7
l2=0.3
wbar=1
r=((1/beta)-1-1e-6 +delta)


## Test code

function f!(x, fvec)
    ps1= wbar + ((kappa*(1-beta*(1-sigma*((1-x[1])/x[1]))))/(beta*((x[1]/(sigma*(1-x[1])))^(gamma/(1-           gamma)))*(1/(2-x[1]))))
    ps2= wbar + ((kappa*(1-beta*(1-sigma*((1-x[2])/x[1]))))/(beta*((x[2]/(sigma*(1-x[2])))^(gamma/(1-gamma)))*(1/(2-x[2]))))

    prod1=prod[1]
    prod2=prod[50]
    y1=(1-x[1])*l1
    y2=(1-x[2])*l2
    M=(((prod1*y1)^((psi-1)/psi))+((prod2*y2)^((psi-1)/psi)))^(psi/(psi-1))
    K=((r/eta)^(1/(eta-1)))*M



    pd1=(1-eta)*(K^eta)*(M^(-eta))*((((prod1*y1)^((psi-1)/psi))+((prod2*y2)^((psi-1)/psi)))^(1/(psi-1)))*    ((prod1*y1)^(-1/psi))*prod1
    pd2=(1-eta)*(K^eta)*(M^(-eta))*((((prod1*y1)^((psi-1)/psi))+((prod2*y2)^((psi-1)/psi)))^(1/(psi-1)))*((prod2*y2)^(-1/psi))*prod2

    fvec[1]=pd1-ps1
    fvec[2]=pd2-ps2
end

mcpsolve(f!,[0.0,0.0],[1.0,1.0], [ 0.3; 0.3])

davidzentlermunro avatar Jul 01 '15 10:07 davidzentlermunro

Getting this also.

From what I can see here it appears that with either the smooth or min max reformulation of the problem the objective is evaluated at the proposed solution before bounds on the proposed solution are enforced.

@sebastien-villemot is this part of the algorithm? is there a way to reformulate the problem such that constraints are satisfied before evaluating the objective?

sglyon avatar Sep 11 '15 13:09 sglyon

Is this still an issue? I am also getting values outside of my defined bounds. I can provide more details if necessary.

jcreinhold avatar Mar 19 '19 02:03 jcreinhold

Is this still an issue? I am also getting values outside of my defined bounds. I can provide more details if necessary.

I don't think anything was done to remedy this, so this is probably still an issue. Is your problem relatively simple?

pkofod avatar Mar 19 '19 08:03 pkofod

Do note the discussion here though: https://github.com/JuliaNLSolvers/NLsolve.jl/issues/100 . The transformation begin done only penalizes deviations from the bounds. The underlying algorithm still assumes that you can evaluate F outside of the bounds. It's not a non-linear solver with constraints as such. I'm not sure what would happen if we just projected the state back into the box, but we may have to be a bit careful if we do.

pkofod avatar Mar 19 '19 08:03 pkofod

The problem isn't too difficult. I am boarding a plane now, but I can post it later. This isn't a huge issue; however, the documentation could be a little clearer about the lack of bound enforcement, for the unenlightened like me. I posted my problem in the forums since it appears better addressed there. Thanks for the package!

jcreinhold avatar Mar 19 '19 10:03 jcreinhold

I'm still having problems with this!

davidzentlermunro avatar Apr 18 '19 17:04 davidzentlermunro