Conditional Constraint Construction Based on Input Data Presence/Absence of Set Combinations in Julia JUMP

I’d write your model quite differently.

Let’s ignore the country/regionn/area/sceario/year/season/time sets for now, because they aren’t relevant to the input/output balance.

You could do something like:

df_carrier = DataFrames.stack(
    CarrierMix, 
    DataFrames.Not([:Direction, :Fuel]);
    variable_name = :Technology,
)
filter!(row -> row.value > 0, df_carrier)
sort!(df_carrier, [:Direction, :Technology])
model = Model()
@variable(model, x[Technology] >= 0)  # Quantity of each technology
@variable(model, inputs[Fuel])  # Quantity of each input fuel
@variable(model, outputs[Fuel])  # Quantity of each output fuel
for (index, df) in pairs(DataFrames.groupby(df_carrier, [:Direction, :Fuel]))
    if index.Direction == "Input"
        @constraint(
            model,
            inputs[index.Fuel] == sum(row.value * x[row.Technology] for row in eachrow(df))
        )
    else
        @constraint(
            model,
            outputs[index.Fuel] == sum(row.value * x[row.Technology] for row in eachrow(df))
        )
    end
end

julia> print(model)
Feasibility
Subject to
 -51 x[Electrolysis] - x[Electric_Storage] - 3 x[Hydrogen_Storage] - 0.17 x[Methanol_Synthesis] - 0.32 x[Ammonia_Synthesis] + inputs[Electricity] = 0
 -33.36111111111111 x[Hydrogen_Storage] - 6.672222222222221 x[Methanol_Synthesis] - 6.08 x[Ammonia_Synthesis] + inputs[Hydrogen] = 0
 -0.84 x[Ammonia_Synthesis] + inputs[Nitrogen] = 0
 -9 x[Electrolysis] + inputs[Water] = 0
 -1.46 x[Methanol_Synthesis] + inputs[CO2] = 0
 -x[Solar_PV] + inputs[Sole] = 0
 -x[Wind_Turbine] + inputs[Wind] = 0
 -5.249999999999999 x[Ammonia_Synthesis] + outputs[Ammonia] = 0
 -7.46 x[Electrolysis] - 2.07 x[Methanol_Synthesis] - 0.26 x[Ammonia_Synthesis] + outputs[Heat] = 0
 -0.95 x[Electric_Storage] - 0.95 x[Solar_PV] - 0.95 x[Wind_Turbine] + outputs[Electricity] = 0
 -33.36111111111111 x[Electrolysis] + outputs[Hydrogen] = 0
 -30.358611111111106 x[Hydrogen_Storage] + outputs[HydrogenCom] = 0
 -5.5361111111111105 x[Methanol_Synthesis] + outputs[Methanol] = 0
 x[Electrolysis] ≥ 0
 x[Electric_Storage] ≥ 0
 x[Hydrogen_Storage] ≥ 0
 x[Methanol_Synthesis] ≥ 0
 x[Solar_PV] ≥ 0
 x[Wind_Turbine] ≥ 0
 x[Ammonia_Synthesis] ≥ 0

My key suggestion would be to stop thinking of your data as 2-d tables with some values of 0 (or missing), and normalize it into some appropriate tabular form. For example, my df_carrier looks like:

julia> df_carrier
23×4 DataFrame
 Row │ Direction  Fuel         Technology          value   
     │ Any        Any          String              Any     
─────┼─────────────────────────────────────────────────────
   1 │ Input      Electricity  Ammonia_Synthesis   0.32
   2 │ Input      Hydrogen     Ammonia_Synthesis   6.08
   3 │ Input      Nitrogen     Ammonia_Synthesis   0.84
   4 │ Input      Electricity  Electric_Storage    1.0
   5 │ Input      Electricity  Electrolysis        51.0
   6 │ Input      Water        Electrolysis        9.0
   7 │ Input      Electricity  Hydrogen_Storage    3.0
   8 │ Input      Hydrogen     Hydrogen_Storage    33.3611
   9 │ Input      Electricity  Methanol_Synthesis  0.17
  10 │ Input      Hydrogen     Methanol_Synthesis  6.67222
  11 │ Input      CO2          Methanol_Synthesis  1.46
  12 │ Input      Sole         Solar_PV            1.0
  13 │ Input      Wind         Wind_Turbine        1.0
  14 │ Output     Ammonia      Ammonia_Synthesis   5.25
  15 │ Output     Heat         Ammonia_Synthesis   0.26
  16 │ Output     Electricity  Electric_Storage    0.95
  17 │ Output     Hydrogen     Electrolysis        33.3611
  18 │ Output     Heat         Electrolysis        7.46
  19 │ Output     HydrogenCom  Hydrogen_Storage    30.3586
  20 │ Output     Methanol     Methanol_Synthesis  5.53611
  21 │ Output     Heat         Methanol_Synthesis  2.07
  22 │ Output     Electricity  Solar_PV            0.95
  23 │ Output     Electricity  Wind_Turbine        0.95

this is a lot easier to work with than the matrix.