diff --git a/.github/workflows/test_and_lint.yml b/.github/workflows/test_and_lint.yml
index 655425d..5019dda 100644
--- a/.github/workflows/test_and_lint.yml
+++ b/.github/workflows/test_and_lint.yml
@@ -17,7 +17,7 @@ jobs:
     - name: Install Dependencies
       shell: bash -el {0}
       run: |
-        conda install ipopt
+        conda install ipopt pip
         python -m pip install --upgrade pip
         pip install .[test]
     - name: Test with pytest
diff --git a/eeco/tests/test_emissions.py b/eeco/tests/test_emissions.py
index e378f96..f04c9e2 100644
--- a/eeco/tests/test_emissions.py
+++ b/eeco/tests/test_emissions.py
@@ -113,8 +113,8 @@ def test_calculate_grid_emissions_cvx(
         emissions_units=emissions_units,
     )
     prob = cp.Problem(cp.Minimize(result), constraints)
-    prob.solve()
-    assert result.value == expected
+    prob.solve(solver=cp.CLARABEL)
+    assert result.value == pytest.approx(expected)
     assert model is None
 
 
diff --git a/eeco/tests/test_utils.py b/eeco/tests/test_utils.py
index c629fbb..7fb76c9 100644
--- a/eeco/tests/test_utils.py
+++ b/eeco/tests/test_utils.py
@@ -53,41 +53,46 @@ def gas_constraint(m, t):
 
 @pytest.mark.skipif(skip_all_tests, reason="Exclude all tests")
 @pytest.mark.parametrize(
-    "consumption_data, varstr1, varstr2, expected",
+    "consumption_data, varstr1, varstr2, time_set, expected",
     [
         (
             {"electric": np.ones(96) * 100, "gas": np.ones(96)},
             "electric",
             "gas",
+            None,
             np.ones(96) * 100,
         ),
+        # non-range (float seconds) index set: the numpy charge_array must be
+        # indexed by position and the Pyomo var by member. Regression for
+        # IndexError under numpy >= 1.24 (float index into a numpy array).
+        (
+            {"electric": np.array([100.0, 200.0, 300.0, 400.0])},
+            "electric",
+            np.array([0.1, 0.2, 0.3, 0.4]),
+            [0.0, 900.0, 1800.0, 2700.0],
+            np.array([10.0, 40.0, 90.0, 160.0]),
+        ),
     ],
 )
-def test_multiply_pyo(consumption_data, varstr1, varstr2, expected):
+def test_multiply_pyo(consumption_data, varstr1, varstr2, time_set, expected):
     model = pyo.ConcreteModel()
-    model.T = len(consumption_data["electric"])
-    model.t = range(model.T)
-    pyo_vars = {}
+    model.T = len(consumption_data[varstr1])
+    model.t = range(model.T) if time_set is None else time_set
+    pos = {t: i for i, t in enumerate(model.t)}
     for key, val in consumption_data.items():
-        var = pyo.Var(model.t, initialize=np.zeros(len(val)), bounds=(0, None))
+        var = pyo.Var(model.t, bounds=(None, None))
         model.add_component(key, var)
-        pyo_vars[key] = var
-
-    @model.Constraint(model.t)
-    def electric_constraint(m, t):
-        return consumption_data["electric"][t] == m.electric[t]
-
-    @model.Constraint(model.t)
-    def gas_constraint(m, t):
-        return consumption_data["gas"][t] == m.gas[t]
+        for t in model.t:
+            var[t].fix(float(val[pos[t]]))
 
     var1 = getattr(model, varstr1)
-    var2 = getattr(model, varstr2)
+    # second operand: a Pyomo var (by name) or a numpy charge_array passed directly
+    var2 = getattr(model, varstr2) if isinstance(varstr2, str) else varstr2
     result, model = ut.multiply(var1, var2, model=model, varstr="test")
     model.objective = pyo.Objective(expr=0)
     solver = pyo.SolverFactory("ipopt")
     solver.solve(model)
-    assert np.allclose([pyo.value(result[i]) for i in range(len(result))], expected)
+    assert np.allclose([pyo.value(result[t]) for t in model.t], expected)
     assert model is not None
 
 
diff --git a/eeco/utils.py b/eeco/utils.py
index b9c0182..49d8c3e 100644
--- a/eeco/utils.py
+++ b/eeco/utils.py
@@ -395,6 +395,13 @@ def multiply(
     ) or isinstance(
         expression2, (SumExpression, IndexedExpression, pyo.Param, pyo.Var)
     ):
+        pos = {t: i for i, t in enumerate(model.t)}
+
+        def at(expr, t):
+            # numpy operands are positional (length T); Pyomo operands are
+            # indexed by time-set member. Map the member to its position.
+            return expr[pos[t]] if isinstance(expr, np.ndarray) else expr[t]
+
         if (not isinstance(expression1, (int, float))) and (len(expression1) > 1):
             if (not isinstance(expression2, (int, float))) and (len(expression2) > 1):
                 # TODO: replace model.t with better way to get dimensions
@@ -402,7 +409,7 @@ def multiply(
                 var = model.find_component(varstr)
 
                 def const_rule(model, t):
-                    return var[t] == expression1[t] * expression2[t]
+                    return var[t] == at(expression1, t) * at(expression2, t)
 
                 constraint = pyo.Constraint(model.t, rule=const_rule)
                 model.add_component(varstr + "_constraint", constraint)
@@ -412,7 +419,7 @@ def const_rule(model, t):
                 var = model.find_component(varstr)
 
                 def const_rule(model, t):
-                    return var[t] == expression1[t] * expression2
+                    return var[t] == at(expression1, t) * expression2
 
                 constraint = pyo.Constraint(model.t, rule=const_rule)
                 model.add_component(varstr + "_constraint", constraint)
@@ -422,7 +429,7 @@ def const_rule(model, t):
             var = model.find_component(varstr)
 
             def const_rule(model, t):
-                return var[t] == expression1 * expression2[t]
+                return var[t] == expression1 * at(expression2, t)
 
             constraint = pyo.Constraint(model.t, rule=const_rule)
             model.add_component(varstr + "_constraint", constraint)