I am currently trying to model a network flow topology optimization in OpenMDAO. One of the design variables is the mass flow from a source into the network. I use mass conservation at every network node to determine the flows in the network. The mass conservation is modeled with an implicit component.
Visulaization of a possible graph structure
The added picture shows an exemplary graph with the source in orange, the targets in blue, and junctions in grey. Overall, there are 70 nodes connected by 72 edges. The mass flow of the source and targets serves as an input to the implicit component, and the 72 edges' mass flows should be determined with 70 residuals.
Below is a minimal example of an implicit component:
import numpy as np
import openmdao.api as om
class MassBalance(om.ImplicitComponent):
def setup(self):
# Declare input variables
self.add_input('m_p', val=1)
self.add_input('m_c', val=np.ones(2))
# Declare output variables
self.add_output('m_i', val=np.ones(5))
# Define Residuals
self.add_residual('m_nodal_balance', shape=(8, 1), desc="Residual for nodal mass balance")
def setup_partials(self):
self.declare_partials(of='m_nodal_balance', wrt=['m_i', 'm_p', 'm_c'])
def apply_nonlinear(self, inputs, outputs, residuals):
m_p = inputs['m_p']
m_c = inputs['m_c']
m_i = outputs['m_i']
residuals['m_nodal_balance'] = a_i.dot(m_i) + a_p.dot(m_p) + a_c.dot(m_c)
def linearize(self, inputs, outputs, partials):
# Declare Options
s = self.options['s']
# Define partials -> mass conservation in every node
partials['m_nodal_balance', 'm_i'] = a_i
partials['m_nodal_balance', 'm_p'] = a_p
partials['m_nodal_balance', 'm_c'] = a_c
a_i = np.array([
[1, 1, 0, 0, 0],
[-1, 0, -1, -1, 0],
[0, 0, 1, 0, -1],
[0, -1, 0, 1, 1]
])
a_p = np.array([[-1], [0], [0], [0]])
a_c = np.array([[0, 0], [1, 0], [0, 1], [0, 0]])
p = om.Problem()
p.model.add_subsystem('MassBalance', MassBalance(), promotes=['*'])
p.setup()
which gives the following error: The number of residuals (8) doesn't match number of outputs (5). If any residuals are added using 'add_residuals', their total size must match the total size of the outputs.
Following the tutorials, I first tried to model the mass conservation with an explicit component and a BalanceComp, but we then got the following error: “If you need a BalanceComp in OpenMDAO with different numbers of state variables and residuals, you'll have to create a custom component.”
If I model the mass conservation as a standalone implicit component, it is unclear to me how we can define the outputs as the residual vector must have the same length as the output vector (even with add_residuals).
I could augment the number of residuals with 0 (or the number of outputs, depending on the network configuration) until the number of residuals matches the number of the output, but that would unnecessarily make the problem bigger. Is there a possibility of modeling an implicit component where the size of the residuals does not match the number of state variables in OpenMDAO, and are there examples in the tutorial that I missed?
I think the core source of confusion is mathematical, and not programtic. The lengths of the variables you've defined in your example don't match up with how the math says they need to be.
To have a well defined system, that can be solved using a linear or nonlinear solver, you must have the same number of state variables as residuals. This likely sounds obvious when states abstractly, but its easy to sometimes get it confused when you write out the equations. You're example has 5 outputs (the nodal mass flows I think) and 3 inputs (split between two variables). Then you define 8 residuals, instead of 5.
However, when I check the dimensions of your actual residual equation I get a result of size 4. So none of the numbers add up to what they should.
Instead I'll make a few guesses as to what your intent was and try to reason from there.
you say " 72 edges' mass flows should be determined with 70 residuals", which violates the rule of an equal number of knowns and unknowns. We have to somehow reduce the problem to only 70 unknown mass-flows. I think 72 comes from the addition of two extra nodes, one as a source and one as a sink. However, the source and sink are not unknowns. You have known values there. So those two nodes are modeled with implicit functions (i.e. no residuals are associated with them), but explicit ones.
This problem is very similar to a node-voltage analysis in electrical circuits, which happens to have an existing OpenMDAO example. In the electrical analysis each node has a voltage that is being solved for (except for the source/sink nodes representing the battery an the ground). You'd need some kind of similar potential type function on the nodes here, which would likely be pressure.
So if there are 70 unknown-nodes, there would be 70 pressure state variables. Then we need to structure a residual. We basically want the net massflow across a node to be 0 at steady state. So 0 = mass_in + mass_out (where I've assumed that either mass_in or mass_out is a negative value).
So conceptually, each edge needs to have its own mass-flow rate. That would be need to be computed based on the relative pressures across the edge (i.e. the state variables at each node).
It looks like your a_i matrix is meant to map to the edges and the shape of that matrix would depend on how many edges you had. So for n_nodes and m_edges, you'd have something of the size (m_edges, n_nodes), which you could multiply by the state vector to get the mass flow across each edge. That would then be a vector of length m_edge, which you could then multiply by another matrix of size (n_nodes, m_deges) that would handle the summation of edge-mass flow for each node.
That way you'd end up with a residual vector of equal length to your state vector. The only tricky bit there is how you handle the boundary conditions of the source and sink nodes. There are different approaches here. Some people like to collapse them into the core matricies themselves. Others prefer to add rows/columns with extra "fake" state variables that get set from the inputs, but which get peeled off before being passed to the residual vector (these two forms are called weak and strong boundary conditions).
So before we can structure an Component for this, we need to get a better definition of the equations you want to solve. If you can express a well defined (i.e. number of states match number of residuals) system of equations, then I can better help you map that to a component. But I believe you need to start with the math first.