Julia problem when using nlsolve -- "MethodError: no method matching"

Question

I am trying to use nlsolve to solve a system of nonlinear equations. I am running the code and function below and receive the error copied below. The output should be a 1x4 matrix. I have put in bold the parts that I think are causing trouble. Thank you very much in advance for any insights.

# Description: this file calibrates trade costs from each country to US
using MAT, Statistics, DataFrames, LinearAlgebra, NLsolve
include("/PATH/sys_tHE_AB.jl")

# parameters
N = 10 # number of possible affiliate host countries
USA = N + 1
ρ = 0.05 # rate of time preference
δ = 0.05 # exogenous death rate (here also endogenous exit)
η = 5.0 # elasticity of substitution
H = η^(-η)*(η-1)^(η-1)
P = ones(1,N+1) # price indexes
n = 1000
T = 31 # number of periods
s = 100
b = 1 # lower bound of the (Pareto) productivity distribution
θ = 3 # shape parameter of the (Pareto) productivity distribution 

# reading in parameters from Matlab files
parameters = matread("parameters.mat")

# creating a structure of the parameters
struct Parameter 
    μ::Matrix{Float64}
    σ::Matrix{Float64}
    ρ::Float64
    δ::Float64
    η::Float64
end
param=Parameter(matread("brownian.mat")["mu_USj"], matread("brownian.mat")["sigma_USj"], ρ, δ, η)
μ_USj=matread("brownian.mat")["mu_USj"]
σ_USj=matread("brownian.mat")["sigma_USj"]
β = sqrt.(2*(ρ + δ)./(σ_USj.^2) .+ (μ_USj./(σ_USj.^2) .- 1/2).^2) .- (μ_USj./(σ_USj.^2) .- 1/2)
α = -sqrt.(2*(ρ + δ)./(σ_USj.^2) .+ (μ_USj./(σ_USj.^2) .- 1/2).^2) .- (μ_USj./(σ_USj.^2) .- 1/2)
disc = ρ + δ .- μ_USj
disc = disc'

# Preallocation
k_e     = zeros(n,N,N+1)
Bo_0    = zeros(n,N,N+1)
Z_HE_0  = zeros(n,N,N+1)
Z_EH_0  = zeros(n,N,N+1)
init_HE           = zeros(n,N,N+1,4)
thresholds_HE_sol = zeros(n,N,N+1,4)
init_OH           = zeros(n,N,4)
thresholds_OH_sol = zeros(n,N,T,4)
status_exp = zeros(n,N,N+1,T)
status     = zeros(n,N,T)

# Useful magnitude to compute the solution of the systems
k_h = H*(1 ./ repeat(z, 1, N)) .^(1-η) .* repeat((P[:,1:10] .^η) .* Q[:,1:10],n)

# Threshold initial conditions
Z_OH_0   = repeat(β[1:N]' ./ (β[1:N].-1)', n, 1) .*  (k_h.^-1) .* (repeat(fh,n,1) + (ρ.+δ.-repeat(μ_USj[1:10]',n,1)) .* repeat(Fh,n,1))
Z_HO_0   = (repeat(β[1:N]' ./ (β[1:N].-1)', n, 1).* (k_h.^-1) .* repeat(fh,n,1)) 
Bo_0_H= k_h ./ repeat(disc[:,1:N],n,1) ./ repeat(β[1:N]',n,1) ./ Z_OH_0.^(repeat(β[1:N]',n,1).-1)
for c in 1:N    
    for d in 1:N+1
        if d != c
            k_e[:,c,d]    = H*(τ[c,d]*w[c]./z./w[d]).^(1-η).*P[d]^η*Q[d] 
            Z_HE_0[:,c,d] = (β[d]/(β[d].-1).*k_e[:,c,d].^-1 .* (fe[c,d] + (ρ+δ-μ_USj[d])*Fe[c,d]))   
            Z_EH_0[:,c,d] = (β[d]/(β[d].-1).*k_e[:,c,d].^-1 .* fe[c,d])        
            Bo_0[:,c,d]   = k_e[:,c,d]./disc[d]./β[d]./(Z_HE_0[:,c,d].^(β[d].-1))
        end
    end
end

**# SOLVING FOR THE THRESHOLDS AND THE PARAMETERS OF THE VALUE FUNCTIONS**
for j in 1:3
   for c in 1:N    
      for d in 1:N+1
         if d != c
             k_e[:,c,d]    = H*(τ[c,d]*w[c]./z./w[d]).^(1-η).*P[d]^η*Q[d] # 1000 x 1 vector
             Z_HE_0[:,c,d] = (β[d]/(β[d].-1).*k_e[:,c,d].^-1 .* (fe[c,d] + (ρ+δ-μ_USj[d])*Fe[c,d])) # 1000 x 1 vector   
             Z_EH_0[:,c,d] = (β[d]/(β[d].-1).*k_e[:,c,d].^-1 .* fe[c,d]) # 1000 x 1 vector        
             Bo_0[:,c,d]   = k_e[:,c,d]./disc[d]./β[d]./(Z_HE_0[:,c,d].^(β[d].-1)) # 1000 x 1 vector
        end
   end
 end

# creating initial Hessian
        for i in 1:length(z)
            if k != j
                init_HE[i,j,k,:] = [Z_HE_0[i,j,k] Z_EH_0[i,j,k] 0 Bo_0[i,j,k]] # 1 x 4 vector
                    if i > 20 && j != 3 && isnan(sum(thresholds_HE_sol[i-1,j,k,:])) == 0
                        init_HE[i,j,k,:] = thresholds_HE_sol[i-1,j,k,:]
                    end
                    if init_HE[i,j,k,3] < 0
                        init_HE[i,j,k,3] = 0.05
                    end
                
                # defining function to be minimized                 
                f_thresholds_HE!(F, Z_AB) = sys_tHE_AB!(F, Z_AB, k_e, fe, Fe, param, j)
                 
                **initial_x = init_HE[i,j,k,:]
                res = nlsolve(f_thresholds_HE!, initial_x, method = :trust_region, xtol=10e-5, ftol=10e-8, iterations=10000, store_trace=true, extended_trace=true)**
            end
        end
end

And the function that it uses is:

# This function computes the system of 4 equations to solve for the export thresholds and parameters B and A
function sys_tHE_AB!(F::AbstractArray, Z_AB::AbstractArray, k_e::Array{Float64, 3}, fe::Matrix{Float64}, Fe::Matrix{Float64}, param::Parameter, j::Int64)  # k_e: (1000, 10, 11)  fe: (11, 11)  Fe:  (11, 11)

    μ = param.μ
    σ = param.σ
    μj = μ[j]
    σj = σ[j]
    ρ = param.ρ
    δ = param.δ

    β =  sqrt(2*(ρ + δ)/(σj^2) + (μj/(σj^2) - 1/2)^2) - (μj/(σj^2) - 1/2)
    α = -sqrt(2*(ρ + δ)/(σj^2) + (μj/(σj^2) - 1/2)^2) - (μj/(σj^2) - 1/2)
    disc = ρ + δ - μj
    
    Y_OE = Z_AB[1]
    Y_EO = Z_AB[2]
    Ae   = Z_AB[3]
    Bo   = Z_AB[4]
    
    F[1] =       Ae*Y_OE^α +        k_e*Y_OE/disc - fe/disc - Fe - Bo*Y_OE^β
    F[2] =       Ae*Y_EO^α +        k_e*Y_EO/disc - fe/disc      - Bo*Y_EO^β
    F[3] =     α*Ae*Y_OE^(α-1) +    k_e/disc    - β*Bo*Y_OE^(β-1)
    F[4] =     α*Ae*Y_EO^(α-1) +    k_e/disc    - β*Bo*Y_EO^(β-1)

    H = [F[1] F[2] F[3] F[4]] 

end

And the error output is:

ERROR: MethodError: no method matching +(::Float64, ::Array{Float64, 3})
For element-wise addition, use broadcasting with dot syntax: scalar .+ array
Closest candidates are:
  +(::Any, ::Any, ::Any, ::Any...) at operators.jl:591
  +(::T, ::T) where T<:Union{Float16, Float32, Float64} at float.jl:383
  +(::Union{Float16, Float32, Float64}, ::BigFloat) at mpfr.jl:414

Answer 1

The problem here has nothing to do with nlsolve. The problem is that within sys_tHE_AB you try to (as the error message says) add a scalar to an array.