Use ModelingToolkit.jl to solve a maximization problem

New to Julia
1

I browsed the website of ModelingToolkit.jl and did not find an example that covers a very simple math problem.

Let’s say that I’d like to maximize y by choosing a value of x:

y = a - b x^2

where a = c + dx - e.

The math is simple. The solution is x^* = \frac{d}{2b} where y is maximized.

I understand that using SymPy.jl is probably the more common approach. However, I’d like to explore ModelingToolkit.jl for its potentially better speed. I’d appreciate if someone could provide an example code to the simple optimization problem above.

2

For this problem, there isn’t much difference in the approach between Symbolics and Sympy, as the derivative is a linear function in the unknown, and Symbolics.solve_for (used below handles those with ease.)

julia> using Symbolics

julia> @variables a b c x
4-element Vector{Num}:
 a
 b
 c
 x

julia> ex = a*x^2 + b*x + c
c + a*(x^2) + b*x

julia> ex′ = Symbolics.derivative(ex, x)
b + 2a*x


julia> Symbolics.solve_for(ex′ ~ 0, x)
(-1//2)*b*(a^-1)

Whereas, you might see:

julia> using SymPy

julia> @syms a b c x
(a, b, c, x)

julia> ex = a*x^2 + b*x + c
   2          
a⋅x  + b⋅x + c

julia> ex′ = diff(ex, x)
2⋅a⋅x + b

julia> solve(ex′, x)
1-element Vector{Sym}:
 -b/(2*a)

or using solveset:

julia> u = solveset(ex′, x)
⎧-b ⎫
⎨───⎬
⎩2⋅a⎭

julia> elements(u)
1-element Vector{Sym}:
 -b/(2*a)

I’m not aware of a direct optimizer in either, but would be happy to pointed to one.

2 Likes
3

That’s the symbolic way, and then symbolic-numeric:

using ModelingToolkit, GalacticOptim

@variables x y
@parameters b c d e
a = c + d*x - e
loss = a - b*x^2
sys = OptimizationSystem(-loss,[x],[b,c,d,e])

u0 = [
    x=>0.0
]
p = [
    b => 1.2
    c => 0.2
    d => 0.3
    e => 0.5
]

prob = OptimizationProblem(sys,u0,p,grad=true)
sol = solve(prob,BFGS())
sol.u # [0.12499999999999999]
3 Likes
4

Thank you very much! At the moment I focus more on the symbolic math. The plan is to run numeric simulation later on, probably in a few weeks (or a few months or years, if I turn out to progress very slowly). Your example codes will be very helpful once I started numerical simulation, I think!

5

Thank you very much!

I see that your codes, there was one equation, ex = a*x^2 + b*x + c, for example.

In the example of y=a−bx^2, where a=c+dx−e, I wonder if it is possible to input two equations in Julia, and ask Julia to solve it?

Of course I could substitute a into y myself, and feed the resulted equation to Julia, using your example codes. My interest is to ask Julia to solve more complex problems where there are more equations involved. That’s why I’m interested in learning how I can feed more than one equations to Julia, and ask Julia to solve them.

6

Sure, both Symbolics and SymPy allow for vectors of equations. Symbolics only has linear equations for now.

An example that would work for either is:

julia> using SymPy

julia> @syms a[1:2,1:2] b[1:2] x y
(Sym[a₁_₁ a₁_₂; a₂_₁ a₂_₂], Sym[b₁, b₂], x, y)

julia> eqs = [Eq(sum(a[i,:] .* (x,y)), b[i]) for i ∈ 1:2]
2-element Vector{Sym}:
 a₁_₁⋅x + a₁_₂⋅y = b₁
 a₂_₁⋅x + a₂_₂⋅y = b₂

julia> solve(eqs, [x,y])
Dict{Any, Any} with 2 entries:
  y => (a₁_₁*b₂ - a₂_₁*b₁)/(a₁_₁*a₂_₂ - a₁_₂*a₂_₁)
  x => (-a₁_₂*b₂ + a₂_₂*b₁)/(a₁_₁*a₂_₂ - a₁_₂*a₂_₁)

Or

julia> using Symbolics

julia> @variables a[1:2,1:2] b[1:2] x y
4-element Vector{Any}:
  Num[a₁ˏ₁ a₁ˏ₂; a₂ˏ₁ a₂ˏ₂]
  Num[b₁, b₂]
 x
 y

julia> eqs = [sum(a[i,:] .* (x,y)) ~ b[i] for i ∈ 1:2]
2-element Vector{Equation}:
 a₁ˏ₁*x + a₁ˏ₂*y ~ b₁
 a₂ˏ₁*x + a₂ˏ₂*y ~ b₂

julia> Symbolics.solve_for(eqs, [x,y])
2-element Vector{Num}:
 -*(a₁ˏ₁^-1)*(a₁ˏ₂*(a₂ˏ₁*b₁*(a₁ˏ₁^-1) - b₂)*((a₁ˏ₂*a₂ˏ₁*(a₁ˏ₁^-1) - a₂ˏ₂)^-1) - b₁)
                         (a₂ˏ₁*b₁*(a₁ˏ₁^-1) - b₂)*((a₁ˏ₂*a₂ˏ₁*(a₁ˏ₁^-1) - a₂ˏ₂)^-1)
3 Likes
7

Thank you!

In general, I have been more interested in packages that are written in Julia, because I think such packages are going to have more speed advantages. However, does it mean that SymPy is at the moment a more capable package than Symbolics, for example, because Symbolics only has linear equations for now?

A side note is that the results from SymPy are much more readable than from Symbolics.