def METHOD_NAME(
    scenario_name, solar_filname=None, use_LP=False, lam=None,
):
    """
    Args:
        scenario_name (str):
            Name of the scenario to create.
        solar_filename (str):
            File containing the solar data.
        use_LP (bool, optional):
            If True, uses LP. Default is False.
        lam (float):
            Value of the dual variable for the chance constraint.
    """
    if solar_filename is None:
        raise ValueError("kwarg `solar_filename` is required")
    if lam is None:
        raise RuntimeError("kwarg `lam` is required")
    data = getData(solar_filename)
    num_scenarios = data['solar'].shape[0]
    scenario_index = extract_scenario_index(scenario_name)
    if (scenario_index < 0) or (scenario_index >= num_scenarios):
        raise RuntimeError('Provided scenario index is invalid (must lie in '
                           '{0,1,...' + str(num_scenarios-1) + '} inclusive)')
    model = pyo.ConcreteModel()
    T   = range(data['T'])
    Tm1 = range(data['T']-1)
    model.y = pyo.Var(T, within=pyo.NonNegativeReals)
    model.p = pyo.Var(T, bounds=(0, data['cMax']))
    model.q = pyo.Var(T, bounds=(0, data['dMax']))
    model.x = pyo.Var(T, bounds=(data['eMin'], data['eMax']))
    if (use_LP):
        model.z = pyo.Var([0], within=pyo.UnitInterval)
    else:
        model.z = pyo.Var([0], within=pyo.Binary)
    ''' "Flow balance" constraints '''
    def flow_balance_constraint_rule(model, t):
        return model.x[t+1]==model.x[t] + \
            data['eff'] * model.p[t] - (1/data['eff']) * model.q[t]
    model.flow_constr = pyo.Constraint(Tm1, rule=flow_balance_constraint_rule)
    ''' Big-M constraints '''
    def big_M_constraint_rule(model, t):
        return model.y[t] - model.q[t] + model.p[t] \
            <= data['solar'][scenario_index,t] + \
            data['M'][scenario_index,t] * model.z[0] # Why indexed??
    model.big_m_constr = pyo.Constraint(T, rule=big_M_constraint_rule)
    ''' Objective function (must be minimization or PH crashes) '''
    model.obj = pyo.Objective(expr=-pyo.dot_product(data['rev'], model.y)
        + data['char'] * pyo.quicksum(model.p)
        + data['disc'] * pyo.quicksum(model.q) + lam * model.z[0],
        sense=pyo.minimize)
    fscr = lambda model: pyo.dot_product(data['rev'], model.y)
    model.first_stage_cost = pyo.Expression(rule=fscr)
    model._mpisppy_node_list = [
        stree.ScenarioNode(name='ROOT', cond_prob=1., stage=1,
            cost_expression=model.first_stage_cost,
            nonant_list=[model.y], scen_model=model)
    ]
    return model