Problem with fixing a value vs constraining a value

fixtest.jl

using JuMP
using Ipopt

const N = 60; #Number of years to calculate (from 2010 onwards)
const tstep = 5; #Years per Period
const ρ = 0.015; #Initial rate of social time preference per year
const θ₂ = 2.8; #Exponent of control cost function

#vvv Currently twitchy when changing out this scale value vvv
const scale = 5/3.666;

gₐ = [0.079*exp.(-0.03*(i-1)) for i in 1:N];
Etree = [3.3*(0.8)^(i-1) for i in 1:N];
rr = [1/((1+ρ)^(tstep*(i-1))) for i in 1:N];
fₑₓ = [if i <= 19 0.25+(1/18)*0.45*(i-1) else 0.7 end for i in 1:N];
pbacktime = [344.0*(0.975)^(i-1) for i in 1:N];

# --- Adjustments ---

const α = 1.45; #Setting alpha or gamma to 0.0 will cause the fix version to be optimal.
const γₑ = 0.3;
#Setting this to 1 will cause the fix solution to solve optimally too.
θ₁ = [pbacktime[i]/θ₂/1000.0 for i in 1:N];
#θ₁ = [1.0 for i in 1:N];

# Setting L[i] == 1 and removing scaling values in YGROSS, CPC and PERIODU equations makes fix optimal
# ---

L = Array{Float64}(undef, N);
L[1] = 6838.0;
A = Array{Float64}(undef, N);
A[1] = 3.8;

for i in 1:N-1
    L[i+1] = L[i]*(10500.0/L[i])^0.134;
    A[i+1] = A[i]/(1-gₐ[i]);
end

optimizer = with_optimizer(Ipopt.Optimizer, print_level=5, max_iter=99900,print_frequency_iter=250,sb="yes")
model = Model(optimizer);

# Rate limit
μ_ubound = [if i < 30 1.0 else 1.2 end for i in 1:N];

# Variables #
@variable(model, 0.0 <= μ[i=1:N] <= μ_ubound[i]); # Emission control rate GHGs
@variable(model, FORC[1:N]); # Increase in radiative forcing (watts per m2 from 1900)
@variable(model, 0.0 <= Tₐₜ[1:N] <= 40.0); # Increase temperature of atmosphere (degrees C from 1900)
@variable(model, -1.0 <= Tₗₒ[1:N] <= 20.0); # Increase temperatureof lower oceans (degrees C from 1900)
@variable(model, Mₐₜ[1:N] >= 10.0); # Carbon concentration increase in atmosphere (GtC from 1750)
@variable(model, Mᵤₚ[1:N] >= 100.0); # Carbon concentration increase in shallow oceans (GtC from 1750)
@variable(model, Mₗₒ[1:N] >= 1000.0); # Carbon concentration increase in lower oceans (GtC from 1750)
@variable(model, E[1:N]); # Total CO2 emissions (GtCO2 per year)
@variable(model, Eind[1:N]); # Industrial emissions (GtCO2 per year)
@variable(model, C[1:N] >= 2.0); # Consumption (trillions 2005 US dollars per year)
@variable(model, K[1:N] >= 1.0); # Capital stock (trillions 2005 US dollars)
@variable(model, CPC[1:N] >= 0.01); #  Per capita consumption (thousands 2005 USD per year)
@variable(model, I[1:N] >= 0.0); # Investment (trillions 2005 USD per year)
@variable(model, S[1:N]); # Gross savings rate as fraction of gross world product
@variable(model, RI[1:N]); # Real interest rate (per annum)
@variable(model, Y[1:N] >= 0.0); # Gross world product net of abatement and damages (trillions 2005 USD per year)
@variable(model, YGROSS[1:N] >= 0.0); # Gross world product GROSS of abatement and damages (trillions 2005 USD per year)
@variable(model, YNET[1:N]); # Output net of damages equation (trillions 2005 USD per year)
@variable(model, DAMAGES[1:N]); # Damages (trillions 2005 USD per year)
@variable(model, Ω[1:N]); # Damages as fraction of gross output
@variable(model, Λ[1:N]); # Cost of emissions reductions  (trillions 2005 USD per year)
@variable(model, MCABATE[1:N]); # Marginal cost of abatement (2005$ per ton CO2)
@variable(model, CCA[1:N] <= 6000.0); # Cumulative industrial carbon emissions (GTC)
@variable(model, PERIODU[1:N]); # One period utility function
@variable(model, CPRICE[1:N]); # Carbon price (2005$ per ton of CO2)
@variable(model, CEMUTOTPER[1:N]); # Period utility
@variable(model, UTILITY); # Welfare function

# Equations #
# Emissions Equation
@constraint(model, [i=1:N], E[i] == Eind[i] + Etree[i]);
# Industrial Emissions
@constraint(model, [i=1:N], Eind[i] == YGROSS[i] * (1-μ[i]));
# Radiative forcing equation
@NLconstraint(model, [i=1:N], FORC[i] == 3.8 * (log(Mₐₜ[i]/588.0)/log(2)) + fₑₓ[i]);
# Equation for damage fraction
@constraint(model, [i=1:N], Ω[i] == 0.0);
# Damage equation
@constraint(model, [i=1:N], DAMAGES[i] == YGROSS[i]*Ω[i]);
# Cost of exissions reductions equation
@NLconstraint(model, [i=1:N], Λ[i] == YGROSS[i] * θ₁[i] * μ[i]^θ₂);
# Equation for MC abatement
@NLconstraint(model, [i=1:N], MCABATE[i] == pbacktime[i] * μ[i]^(θ₂-1));
# Carbon price equation from abatement
@NLconstraint(model, [i=1:N], CPRICE[i] == pbacktime[i] * μ[i]^(θ₂-1));
# Output gross equation
@NLconstraint(model, [i=1:N], YGROSS[i] == A[i]*(L[i]/1000.0)^(1-γₑ)*K[i]^γₑ);
# Output net of damages equation
@constraint(model, [i=1:N], YNET[i] == YGROSS[i]*(1-Ω[i]));
# Output net equation
@constraint(model, [i=1:N], Y[i] == YNET[i] - Λ[i]);
# Consumption equation
@constraint(model, [i=1:N], C[i] == Y[i] - I[i]);
# Per capita consumption definition
@constraint(model, [i=1:N], CPC[i] == 1000.0 * C[i] / L[i]);
# Savings rate equation
@constraint(model, [i=1:N], I[i] == S[i] * Y[i]);
# Period utility
@constraint(model, [i=1:N], CEMUTOTPER[i] == PERIODU[i] * L[i] * rr[i]);
# Instantaneous utility function equation
@NLconstraint(model, [i=1:N], PERIODU[i] == ((C[i]*1000.0/L[i])^(1-α)-1)/(1-α)-1);

# Equations (offset) #
# Cumulative carbon emissions
@constraint(model, [i=1:N-1], CCA[i+1] == CCA[i] + Eind[i]*scale);
# Atmospheric concentration equation
@constraint(model, [i=1:N-1], Mₐₜ[i+1] == Mₐₜ[i]*0.912 + Mᵤₚ[i]*0.03832888888888889 + E[i]*scale);
# Lower ocean concentration
@constraint(model, [i=1:N-1], Mₗₒ[i+1] == Mₗₒ[i]*0.9996625 + Mᵤₚ[i]*0.0025);
# Shallow ocean concentration
@constraint(model, [i=1:N-1], Mᵤₚ[i+1] == Mₐₜ[i]*0.088 + Mᵤₚ[i]*0.9591711111111112 + Mₗₒ[i]*0.0003375);
# Temperature-climate equation for atmosphere
@constraint(model, [i=1:N-1], Tₐₜ[i+1] == Tₐₜ[i] + 0.098*((FORC[i+1]-(3.8/2.9)*Tₐₜ[i])-(0.088*(Tₐₜ[i]-Tₗₒ[i]))));
# Temperature-climate equation for lower oceans
@constraint(model, [i=1:N-1], Tₗₒ[i+1] == Tₗₒ[i] + 0.025*(Tₐₜ[i]-Tₗₒ[i]));
# Capital balance equation
@NLconstraint(model, [i=1:N-1], K[i+1] <= 0.9^tstep * K[i] + tstep*I[i]);
# Interest rate equation
@NLconstraint(model, [i=1:N-1], RI[i] == (1+ρ)*(CPC[i+1]/CPC[i])^(α/tstep)-1);

# Savings rate for asympotic equilibrium
for i=N-9:N
    JuMP.fix(S[i], 0.2582782; force=true);
end
# Initial conditions
JuMP.fix(CCA[1], 90.0; force=true);
JuMP.fix(K[1], 135.0; force=true);
JuMP.fix(Mₐₜ[1], 830.4; force=true);
JuMP.fix(Mᵤₚ[1], 1527.0; force=true);
JuMP.fix(Mₗₒ[1], 10010.0; force=true);
JuMP.fix(Tₐₜ[1], 0.8; force=true);
JuMP.fix(Tₗₒ[1], 0.0068; force=true);

@constraint(model, UTILITY == tstep * 0.016408662 * sum(CEMUTOTPER[i] for i=1:N) - 3855.106895);

# Objective function
@objective(model, Max, UTILITY);

# -- MODIFICATION --

#Here's the modification I suspect should be right:
for i in 2:N
    JuMP.set_upper_bound(Tₐₜ[i], 2.0);
end

# Here's the solution that works.
#JuMP.unfix(Tₐₜ[1]);
#JuMP.set_lower_bound(Tₐₜ[1], 0.0);
#@constraint(model, Tₐₜ[1] == 0.8);
#for i in 1:N
#   JuMP.set_upper_bound(Tₐₜ[i], 2.0);
#end

optimize!(model);