Any concise way to express simulation equations in Julia (like algebraic modelling language but not necessarily linked to optimisation problems)?

Specific Domains Data
1

I know several Algebraic Modelling Languages in the context of optimisation, e.g. GAMS, Pyomo (Python) or JuMP (Julia). These allow you to express your objective function and your constraints with a syntax very concise and similar to your pen and paper equations.
But… outside the domain of optimisation, when you have a model with hundred of equations like the following pseudocode, which is the simplest way to express and compute these ?

@meq nTrees!(r in reg, sp in species, dc in diameterClass[2-end], y in years) = nTrees_(r, sp, dc, y-1)*(1-mortRate_(r, sp, dc, y-1) - promotionRate_(r, sp, dc, y-1))) +  nTrees_(r, sp, dc-1, y-1) * promotionRate_(r, sp, dc-1, y-1))
....
production!(r in reg, sp in secPr, y in years)   = sum(trValues_(r,pp) * transfCoef_(r,pp)  for pp in primPr)

I assume a matrix formulation of the equations or maybe using for loops/list comprehension/maps are the most efficient solutions, but they become cumbersome when you have many dimensions and you need to operate on a subset of them, or when you have to operates transforms like the summation above…

What is your approach? Is there a library that allows to write the kind of equations above in a concise way ? Is such library needed ? (I did wrote one some time ago based on IndexedTable.jl and have to decide if it is worth to update it to Julia 1.x or if there are better approaches)

2

So far it covers:

  • ODEs
  • SDEs
  • Jump dynamics
  • Nonlinear solving
  • Nonlinear optimization
  • Control
  • Reaction systems
  • PDEs
2 Likes
3

Thank you, but my needs are more on how to manipulate data (likely one of the two commonly used data structures in Julia, IndexedTables or DataFrames) in a concise way…

4

I would use an inner join for the first example and split-apply-combine for the second.
There are also several data manipulation packages like Query.jl that might allow a more concise syntax for the same operations.

5

Thank you, but even using Query.jl, it would take lot of efforts to write every single equation/transformation.

Here is an extract of the model I am playing with, using the (old, not working) MultiDimEquation.jl package I am considering to revive or looking if something similar (and better) exists already. Note that the functions are auto-generated with Foo_(dim1,dim2) a get_like function and Foo!(value,dim1,dim2) a set-like function:

# *** Init market module ***

# Setting da for prim products and sa for sec products for each region from the data at national level
@meq  sa!(r in fr2, tp in secProducts,y2) = sa_("FRA",tp,y2)
@meq  da!(r in fr2, pp in priProducts,y2) = da_("FRA",pp,y2)

# set local prices for transformed products
@meq  pl!(r in fr2, tp in secProducts,y2) = sum(pl_(r,pp,y2)*a_("FRA",pp,y1,tp) for pp in priProducts) + m_(r,tp,y1)

# set local demand for primary products
# [dl!( sum(sl_(r,tp) * a_(r,pp,tp) for tp in secProducts)  ,r, pp) for r in fr2, pp in priProducts ]
@meq dl!(r in fr2, pp in priProducts, y2) = sum(sl_(r,tp,y2) * a_("FRA",pp,y1,tp) for tp in secProducts)

# set total supply (local + abroad[exports]) and total demand (local + abroad[imports])
@meq st!(r in fr2, p in products, y2) = sl_(r, p, y2) + sa_(r, p, y2)
@meq dt!(r in fr2, p in products, y2) = da_(r, p, y2) + dl_(r, p, y2)

# set initial values for the ces function (transf products)
@meq q1!(r in fr2, tp in secProducts,y1) = 1.0 /  (1.0+((dl_(r,tp,y2)/da_(r,tp,y2))^(1/psi_(r,tp,y1)) ) * (pl_(r,tp,y2)/pl_(wld,tp,y2))   )
@meq p1!(r in fr2, tp in secProducts,y1) = 1 - q1_(r, tp,y1)

# composite demand/price (sec products)
@meq dc!(r in fr2, tp in secProducts,y2) = (q1_(r,tp,y1)*da_(r,tp,y2)^( (psi_(r,tp,y1)-1)/psi_(r,tp,y1)  )+p1_(r,tp,y1)*dl_(r,tp,y1)^( (psi_(r,tp,y1)-1)/psi_(r,tp,y1) )  )^(psi_(r,tp,y1)/(psi_(r,tp,y1)-1))
@meq pc!(r in fr2, tp in secProducts,y2) = (da_(r,tp,y2)/dc_(r,tp,y2)) * pl_(wld,tp,y2) + (dl_(r,tp,y2)/dc_(r,tp,y2)) * pl_(r,tp,y2)

# Capacity parameter
@meq k!(r in fr2, tp in secProducts,y2) = k1_("FRA",tp,y1) * sl_(r,tp,y2)

# set initial values for the ces function (primary products)
@meq t1!(r in fr2, pp in priProducts,y1) = 1.0 /  (1.0+((sl_(r,pp,y2)/sa_(r,pp,y2))^(1/eta_(r,pp,y1)) ) * (pl_(r,pp,y2)/pl_(wld,pp,y2))   )
@meq r1!(r in fr2, pp in priProducts,y1) = 1 - t1_(r, pp,y1)

# composite demand/price (prim products)
@meq sc!(r in fr2, pp in priProducts,y2) = (t1_(r,pp,y1)*sa_(r,pp,y2)^( (eta_(r,pp,y1)-1)/eta_(r,pp,y1)  )+r1_(r,pp,y1)*sl_(r,pp,y2)^( (eta_(r,pp,y1)-1)/eta_(r,pp,y1) )  )^(eta_(r,pp,y1)/(eta_(r,pp,y1)-1))
@meq pc!(r in fr2, pp in priProducts, y2) = (sa_(r,pp,y2)/sc_(r,pp,y2)) * pl_(wld,pp,y2) + (sl_(r,pp,y2)/sc_(r,pp,y2)) * pl_(r,pp,y2)
# psi_("LO","hardWRoundW",y1)

# weighted prices
@meq pw!(r in fr2, tp in secProducts,y2) = (dl_(r,tp,y2)*pl_(r,tp,y2) + da_(r,tp,y2)*pl_(wld,tp,y2))/ (dl_(r,tp,y2)+da_(r,tp,y2))
# at this point there is an error, as some outputs are NaN!!
@meq pw!(r in fr2, pp in priProducts,y2) = (sl_(r,pp,y2)*pl_(r,pp,y2) + sa_(r,pp,y2)*pl_(wld,pp,y2))/(sl_(r,pp,y2)+sa_(r,pp,y2))

# setting regional trade to zero and positive transport costs between same region
@meq rt!(r in fr2, p in products, y2, r2 in fr2) = 0.0
[ct!(0.1, rFrom, p, y1, rTo ) for rFrom in fr2, p in products, rTo in fr2 if rFrom == rTo ]

# available inventories
@meq in!(r in fr2, pp in priProducts, y2) = sum(app_(pp,ft,dc)* vol_(r,ft,dc,y1) for ft in fTypes, dc in invDClasses)
# [println("$r \t $p \t $(in_(r,p,y2))") for r in fr2, p in priProducts if r == "AQ"]

# *** Forest dynamic module ***

# Harvesting rate
@meq hr!(r in fr2, ft in fTypes, dc in invDClasses, y2) = sum(app_(pp,ft,dc) * st_(r,pp,y2) / in_(r, pp, y2) for pp in priProducts)
@meq hV!(r in fr2, ft in fTypes, dc in invDClasses, y2) = hr_(r, ft, dc, y2) * vol_(r, ft, dc, y1)

# volumes.. first diameter class and then subsequent diameter classes..
@meq (vol!(r in fr2, ft in fTypes, invDClasses[1], y2) =
           vol_(r, ft, invDClasses[1], y1) * (1 - mortCoef_(r, ft, invDClasses[1], y1) - (1/tp_(r, ft, invDClasses[1], y1)) -hr_(r, ft, invDClasses[1], y2))
        +  vReg_(r, ft, invDClasses[1], y1))

@meq  (vol!(r in fr2, ft in fTypes, dc in invDClasses[2:end], y2) =
           vol_(r, ft, dc, y1) * (1-mortCoef_(r, ft, dc, y1) - 1/tp_(r, ft, dc, y1)-hr_(r, ft, dc, y2))
        +  vol_(r, ft, dc-1, y1) * (1/tp_(r, ft, dc-1, y1)) * betaCoef_(r, ft, dc-1, y1))