Constrained Optimization of a "taxation" curve not working Optimization.jl / Ipopt

Specific Domains Optimization (Mathematical)
,
1

I was having some discussions online on Mastodon, and advocating that we could simplify our tax code to a flat tax + UBI and that this is kind of “near optimal” in some sense.

I wrote up a notebook about it… GitHub - dlakelan/UBINotebook: A notebook and articles on the advantages of Universal Basic Income

I started out using IPNewton from Optim.jl and got some results, but was having some problems and I decided to try Ipopt through OptimizationMOI.jl

What I’m doing is constructing an Approxfun Fun object from a list of coefficients, and then evaluating the quality of this function to get a “fitness” and trying to optimize this fitness subject to two constraints (and then following up later with maybe some more constraints)

Right now it’s stopping with a ridiculous solution that has vastly violated constraints. You can check out the github repo to see what’s up.

Any hints on the best solver to use to get general equality and inequality constraints?

EDIT: the section where the optimization is supposed to happen is like this:

using DataFrames, DataFramesMeta, StatsPlots, Distributions, ApproxFun, Optimization, OptimizationBase, 
    OptimizationMOI, Ipopt,
    Printf, Interact, Conda

const nhous = 125736353.0
const currevenue = 3.8 # current revenue from income and payroll tax in trillions

incDist = MixtureModel([Uniform(0,28e3),Uniform(28e3,55e3),Uniform(55e3,89.7e3),Uniform(89.7e3,149e3), Truncated(Exponential(269356),115e3,2e6)])

incsamp = rand(incDist,10000)


function coefstofun(u,p)
    spc = Chebyshev(0.0..2_000_000.0)
    fsqt = Fun(spc,u)
#    df = fsqt*fsqt
    df = fsqt
    fI = Integral() * df
    ffI = fI - fI(0.0) + p.min
    return ffI,df
end

function fitness(u,p)
    fI,df = coefstofun(u,p)
    incomes = p.incs
    
    meaninc = mean(fI(i) for i in incomes)
    meanmarg = mean(df(i) for i in incomes)
    fitness = -(meaninc/p.gdppc + p.k * meanmarg)
end

function revconstr(constr, u,p) ## we need a certain revenue
    target = p.revneeded # in trillions
    fI,df = coefstofun(u,p)
    revactual = mean(i - fI(i) for i in p.incs ) * (nhous / 1e12) 
    constr[1] = revactual - target 
    return constr[1]
end

function revderivconstr(constr,u,p)
    target = p.revneeded
    incs = p.incs
    fI,df = coefstofun(u,p)
    revactual = mean(i-fI(i) for i in incs) * (nhous/1e12)
    constr[1] = revactual - target
    mind = minimum(df(i) for i in incs)
    constr[2] = mind
    return constr
end

let #u0 = [0.44678816662151877, -0.08788321820882686, 0.17815799706452318, -0.14584610993899744, 0.07437859490773412, -0.12028356099198997, -0.27147771772117185]
    u0 = rand(Normal(0.0,0.0001),7) .+ [.6; zeros(6)]

@manipulate for min=0.01:.01:.25, k=0.0:0.01:1.0, maxtax = 0.2:0.01:0.95
    gdpc = mean(incsamp)
    parms = (k=k,min=min*gdpc,incs=incsamp,revneeded=currevenue,gdppc=gdpc)
    tax = (currevenue + min*gdpc*nhous/1e12) / (gdpc * nhous/1e12)
    u0 = [1.0- 1.03*tax,0.01,0.01,0.0,0.0,0.0,0.0]
    earn = collect(0.0:1000.0:1000000.0)
    p = plot(earn,earn,xlim=(0.0,400e3),ylim=(0.0,400e3),label="y = x visual reference",xlab="Earned Income (Dollars)", ylab="Take Home Income (Dollars)")
    p = plot!(earn,[min*gdpc + e*(1.0-tax) for e in earn],label="flat tax+UBI")
    opprob = OptimizationProblem(OptimizationFunction(fitness,Optimization.AutoForwardDiff(),cons=revderivconstr),
            u0,parms, 
            lcons=[0.0,1.0-maxtax],ucons=[0.0,3.0])
    sol = solve(opprob,Ipopt.Optimizer(); max_wall_time=35.0, print_level=5,output_file="optim.out")
#    println(sol)
    revcon=[0.0, 0.0]
    revderivconstr(revcon,sol.u,parms)
    #u0 = sol.u .+ rand(Normal(0.0,.04),7)
#    println("Revenue constraint: $(revcon[1])")
#    println("Deriv constraint: $(revcon[2])")
    fI,df = coefstofun(sol.u,parms)
    p = plot!(earn,[fI(e) for e in earn],label="Optimal Tax Curve")
    
end

end
2

Can you include the various using Ipopt etc so that the example is runnable?

Edit: I guess there’s the notebook instead

3

Just based on the output in the notebook,

image
image2334×876 239 KB

The solution violates the constraints because Ipopt did not find a solution. It hit the time limit instead.

4

I added those things, but it’s probably easier to check out the github? but let me know if you can make it work with what I put.

5

Can you simplify as much as possible? I likely don’t need to install Conda, and I don’t need the Interact slider. I just want to be able to copy-paste one block of code and hit the problematic solve.

6

Sure. I will have to come back to this in an hour or so as family dinner is begin served. I’ll try to strip out the part involving the interactive bits etc. The goal of the thing was to let people run the notebook, but I agree to debug the optimization doesn’t need all that.

I will say IPNewton was able to solve this problem in 5 seconds while Ipopt is not liking it even after 35 seconds… when I let it go as long as it wanted it gave a ridiculous and infeasible solution.

7

Ok, I’m pushing a branch “ipoptdebug” which has ipoptdebug.jl a script you can just include… I’m waiting for it to precompile all the stuff so I can test it… it’s taking a while

8 
using ApproxFun
import Distributions
import Ipopt
import Optimization
import OptimizationMOI
import Statistics

function coefs_to_fun(u, p)
    spc = ApproxFun.Chebyshev(0.0..2_000_000.0)
    fsqt = ApproxFun.Fun(spc, u)
    fI = ApproxFun.Integral() * fsqt
    ffI = fI - fI(0.0) + p.min
    return ffI, fsqt
end

function fitness(u,p)
    fI, df = coefs_to_fun(u, p)
    meaninc = Statistics.mean(fI(i) for i in p.incs)
    meanmarg = Statistics.mean(df(i) for i in p.incs)
    return -(meaninc / p.gdppc + 0.5 * meanmarg)
end

function revderivconstr(constr, u, p)
    fI, df = coefs_to_fun(u, p)
    revactual = Statistics.mean(i - fI(i) for i in p.incs) * (nhous / 1e12)
    constr .= [revactual - p.revneeded, minimum(df(i) for i in p.incs)]
    return constr
end

nhous = 125736353.0
currevenue = 3.8

incDist = Distributions.MixtureModel([
    Distributions.Uniform(0, 28e3),
    Distributions.Uniform(28e3, 55e3),
    Distributions.Uniform(55e3, 89.7e3),
    Distributions.Uniform(89.7e3,149e3),
    Distributions.Truncated(Distributions.Exponential(269356), 115e3, 2e6),
])
incsamp = rand(incDist, 10_000)
min = 0.1
gdpc = Statistics.mean(incsamp)
parms = (
    min = min * gdpc,
    incs = incsamp,
    revneeded = currevenue,
    gdppc = gdpc,
)
tax = (currevenue + min * gdpc * nhous / 1e12) / (gdpc * nhous / 1e12)
u0 = [1.0 - 1.03 * tax, 0.01, 0.01, 0.0, 0.0, 0.0, 0.0]
opprob = Optimization.OptimizationProblem(
    Optimization.OptimizationFunction(
        fitness,
        Optimization.AutoForwardDiff(),
        cons = revderivconstr,
    ),
    u0,
    parms,
    lcons = [0.0, 0.5],
    ucons = [0.0, 3.0],
)
sol = Optimization.solve(opprob, Ipopt.Optimizer(); max_wall_time = 35.0)
1 Like
9

I’m not sure why, I’d have thought the precompile occurred when I ran the notebook, but now that I’m running from the terminal everything is recompiling at the moment… it’s taking many tens of minutes… sigh

10

Ok, yay I got it all compiled up and it runs. The version I have is similar to yours with less formatting changes. It’s checked in on the ipoptdebug branch at github.

The result I get is:

EXIT: Iterates diverging; problem might be unbounded.
retcode: Failure
u: [6.637830412091017e19, 4.571391715648334e19, -7.873474346436485e19, -4.8963343822844094e20, -1.4371180114549217e19, -2.0462841609132576e20, 2.6619329494434972e20]
Final objective value:     -6.68173093466959e18

Revenue constraint: 1.1565319951125716e6
Deriv constraint: -5.567846983342287e20

The Deriv constraint says the minimum derivative at the sample points should be between 1-maxtax and 10, maxtax is 0.7 so at every sample point from the income distribution the curve should have derivative between 0.3 and 10 instead it has -5.5678e20

The initial conditions might violate the equality constraint (total revenue must be some number, so the difference must be 0.0) but the derivative being between 0.3 and 10 is satisfied by the initial conditions.

Question: do the constraints need to be differentiable? Because one of my constraints is basically that the minimum derivative across the sample points is between 0.3 and 10, and taking the minimum is not really a differentiable thing. So is that screwing with the algorithm?

If the curve has a derivative that is everywhere greater than 0 and starts at a value above 0 then it is an increasing function that’s everywhere positive. Instead the optimizer seems to be finding a function where we make the take-home income of rich people -1e20 dollars then we can give 1e18 dollars to each poor person or something dumb like that… violating the constraint makes the problem diverge.

11

https://coin-or.github.io/Ipopt/index.html#Overview

Says that the constraints should be twice continuously differentiable, so I think my minimum constraint is the issue. I will see about reformulating the constraint. Perhaps I should just make the value of the derivative at each of the chebyshev points be a separate constraint. Or sum up exp(-dy/dx) and make it be less than some number.

12

All right, I tried a constraint in which I calculate the mean of exp(-20*dy/dx) and constrain its upper bound. It works sufficiently well that I can let people tune the peak marginal tax rate by controlling the upper bound, and it converges in a couple seconds so everything is working. It comes down to the need for a twice differentiable constraint for Ipopt.