Design gear train with specific ratio with a set of fixed gears

Specific Domains Optimization (Mathematical)
1

I want to make a gear train out of lego gears that needs a specific gear ratio. So let’s say I want a gear train with a ratio of 45/28 or whatever, and I have gears with [8, 12, 16, 20, 24, 36, 40] teeth. What combination will give me the specified ratio?

It seems like some sort of factorization or knapsack problem, but I’m not sure if there is any algo out there for this particular case.

I’ve defined all the gear ratios as follows

gears = [8, 12, 16, 20, 24, 36, 40]
gearratios = Set()
for g1 = gears
    for g2 = gears
        if g1 > g2
            push!(gearratios, g1/g2)
        end
    end
end

And then I make an array of integers and have a cost function like

cost(x) = abs(prod(gearratios.^x)/ratio-1)

I was thinking, maybe just run some optimizer on it, but I’m not sure which kind would deal well with this type of problem. JuMP wasn’t happy wiith the exponent, and messing around with simulated annealing didn’t give any useful results so far. For that I was using a neighbor function that just increments/decrements a random number in x.

Uh, ideas and suggestions welcome.

1 Like
2

You should probably tell the optimizer that x is vector of Ints. You might be able to keep everything in rational arithmetic if you minimize

abs( prod(gearrations.^x) - ratio )

instead.

3

Yea I tried that the question is more, which optimizer can deal with integer arrays and this cost function?

For example, JuMP gives the following error

@objective(model, Min, abs(prod(gearratios.^x)-ratio))

ERROR: MethodError: no method matching ^(::Rational{Int64}, ::VariableRef)

So it seems you just can’t express this problem in JuMP.

4

After a bit more experimenting, I sorta got simulated annealing to work.

using Optim

M2 = 12.4206012
S2 = 12
N2 = 12.65834751
K1 = 23.93447213
O1 = 25.81933871

periods = [M2, S2, N2, K1, O1]./12

ratios = rationalize.(periods; tol=0.001)

gears = [8, 12, 16, 20, 24, 36, 40]
gearratios = Set()
for g1 = gears
    for g2 = gears
        if g1 > g2
            push!(gearratios, g1//g2)
        end
    end
end

x0 = zeros(length(gearratios))


function neighbor(x::AbstractArray, x_proposal::AbstractArray)
    for idx =1:length(x)
        val = rand((1, 0, -1))
        x_proposal[idx] = x[idx]+val
        # println(x_proposal)
    end
end

for ratio = ratios
    cost(x) = abs(prod(gearratios.^x)-ratio)#+sum(abs.(x))/1e3
    res = optimize(cost, x0, SimulatedAnnealing(;neighbor))
    println()
    println(ratio)
    println(Optim.minimum(res))
    println(Optim.minimizer(res))
end

But of course, by only optimizing for the error it will happily produce solutions that require tons and tons of gears.

29//28
0.012684630641507377
[0.0, 6.0, 0.0, -8.0, 0.0, -10.0, 4.0, 8.0, 0.0, -5.0, -4.0, 1.0, 6.0, -5.0]

1//1
0.0
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]

19//18
0.00169363913017051
[4.0, 0.0, -3.0, 0.0, -1.0, 3.0, -5.0, 1.0, 3.0, -1.0, 5.0, -3.0, -10.0, 8.0]

309//155
0.02916766228594181
[-2.0, -6.0, 2.0, 2.0, -3.0, 0.0, 3.0, -4.0, 1.0, 6.0, 5.0, 1.0, 5.0, 0.0]

71//33
0.012000565780664996
[3.0, -1.0, 5.0, 2.0, 1.0, -1.0, -4.0, -3.0, 0.0, 1.0, -2.0, 2.0, -6.0, -1.0]
5

You could minimize the number of gears with prod(gearrations.^x) = ratio as a constraint.

6

That assumes that the constraint is easy to satisfy, and so far none of my solutions even get the exact result. I think it’s probably easier to constrain the number of gears and then optimize for error.

7

Some ideas (withouth implementing them in Julia). Perhaps the logarithm can help here a bit. Let me rephrase the problem, just to make sure: you want find a set of nonnegative integers x_1, x_2, \ldots, x_n, such that

r_1^{x_1}r_2^{x_2}\ldots r_n^{x_n} = d,

where d is the demanded total ratio, n is the number of the (kinds of) gear ratios available to you and r_i are the individual gear ratios (you can have more than one piece of each kind).

You can apply logarithm function to the constraint equation to get

\ln(r_1^{x_1}r_2^{x_2}\ldots r_n^{x_n}) = \ln d,

which is equivalent to

x_1\ln r_1 + x_2\ln r_2 + \ldots + x_n\ln r_n= \ln d,

which after relabelling is

a_1x_1 + a_2x_2 + \ldots + a_nx_n= b.

Of course, the crucial challenge is that x_i s are nonnegative integers. If only we were guaranteed that a solution exists, we could formulated the following integer optimization problem

minimize x_1 + x_2 + \ldots + x_n
over x_i\in Z_+, i=1,\ldots, n
subject to a_1x_1 + a_2x_2 + \ldots + a_nx_n= b,

which could be rewritten compactly as

minimize \|x\|_1
over x\in Z_+^n
subject to a^Tx= b.

But generally an existence of a solution will hardly be guaranteeable here. Perhaps reformulate then to something like

minimize (a^Tx- b)^2+ \gamma\|x\|_1
over x\in Z_+^n

if at least approximate solution (composed of as few as possible gears, sort of, it is not true minimization, just regularization) is acceptable.

4 Likes
8

I am giving it a try with some Julia code too:

g = [8, 12, 16, 20, 24, 36, 40]                     # Available gears.
r = []                                              # Initialize an array of gear ratios.
for g1 in g                                         # Fill in the array of gear ratios with values.
    for g2 in g
        if g1 > g2                                  # If nonnegativity constraints on the xᵢ variables are to be imposed below, remove the condition; 
            push!(r, g1/g2)                         #  otherwise only gear ratios > 1 can be achieved.   
        end
    end
end

a = log.(r)
n = length(a)                                       # Number of (types of) gear ratios.

d = 45/28                                           # Required gear ratio.
b = log(d)

using JuMP
using LinearAlgebra
using Juniper
using Ipopt

optimizer = Juniper.Optimizer
nl_solver = optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0)

optimizer = Juniper.Optimizer
model = Model(optimizer_with_attributes(optimizer, "nl_solver"=>nl_solver))
@variable(model, x[1:n]>=0,Int)                     # If assumed nonnegative, only gear ratios >=1 can be obtained. 
#@variable(model, x[1:n],Int)                       # But currently fails if the nonnegativity constraint is removed. Dunno why.

γ = 0.001                                           # Regularization coefficient (if small, it does not try to minimize the number of gears ratios)
@NLexpression(model,expr,(sum(a[i]*x[i] for i = 1:n)-b)^2 + γ*sum(abs(x[i]) for i=1:n))
@NLobjective(model, Min, expr)
print(model)

optimize!(model)

It finally (originally I reported troubles here) seems doing what I expected:

julia> println(JuMP.value.(x))
[0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]

julia> println(JuMP.objective_value(model))
0.0020198412370350254

julia> println(JuMP.termination_status(model))

Let’s compare the achieved gear ratio with the demanded one:

julia> r[3]*r[10]
1.5999999999999999

julia> d
1.6071428571428572

Sort of…

Is this satisfactory, @pepijndevos ?