JuliaNLSolvers / NLsolve.jl

Julia solvers for systems of nonlinear equations and mixed complementarity problems
Other
330 stars 66 forks source link

Function straying out of bounds in NLsolve #19

Open davidzentlermunro opened 9 years ago

davidzentlermunro commented 9 years ago

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])
sglyon commented 9 years ago

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?

jcreinhold commented 5 years ago

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

pkofod commented 5 years ago

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 commented 5 years ago

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.

jcreinhold commented 5 years ago

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!

davidzentlermunro commented 5 years ago

I'm still having problems with this!