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

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)

So far it covers:

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

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…

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.

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))
``````