---
name: optiskill
description: >
  Interactive dialogue system for querying, analysing, and manipulating Pyomo
  optimization models using natural language. Supports six analysis modes:
  model description, retrieval, sensitivity (dual values), what-if (parameter
  changes), why-not (forced decisions), and custom Python code execution.
  Trigger on any question about an optimization model: "What does this model do?",
  "What is the value of X?", "How sensitive is cost to demand?", "What if we
  double capacity?", "Why doesn't the model choose route A?", "Add a new
  constraint", "Extract dual values".
allowed-tools: Bash(python* pip*), Read, Write
---

# OptiSkill — Pyomo Optimization Model Analysis

## Goal

Answer natural-language questions about Pyomo optimization models by classifying
each query into one of six analysis modes and executing the appropriate workflow.

---

## Prerequisites — Environment Setup

> Run these steps **once per session**. If already set up, skip to
> **Query Classification**.

### 0a. Install dependencies

```bash
pip install pyomo cloudpickle highspy
```

> Fallback: if `highspy` fails, install GLPK:
> `apt-get install -y glpk-utils` or `conda install -c conda-forge glpk`

### 0b. Write the model source and data

```bash
mkdir -p /tmp/optiskill_workspace
```

<details>
<summary>diet.py — Stigler Diet Problem (click to expand)</summary>

```bash
cat > /tmp/optiskill_workspace/diet.py << 'PYEOF'
import json
from pyomo.environ import *

model = ConcreteModel()
data = globals().get("data", {})
nutrient_requirements = data['nutrient_requirements']
nutritive_values = data['nutritive_values']

model.n = Set(initialize=list(data['nutrient_requirements'].keys()), doc='nutrients')
model.f = Set(initialize=list(data['nutritive_values'].keys()), doc='foods')
model.b = Param(model.n, initialize=nutrient_requirements, mutable=True,
                doc='required daily allowances of nutrients')
initialized_data = {
    (food, nutrient): value
    for food, nutrients in nutritive_values.items()
    for nutrient, value in nutrients.items()
}
model.a = Param(model.f, model.n, initialize=initialized_data, mutable=True,
                doc='nutritive value of foods (per dollar spent)')
model.x = Var(model.f, within=NonNegativeReals, bounds=(0, None), initialize=0,
              doc='dollars of food f to be purchased daily (dollars)')
model.cost = Var(within=NonNegativeReals, bounds=(0, None), initialize=0,
                 doc='total food bill (dollars)')

def nutrient_balance_rule(model, n):
    return sum(model.a[f, n] * model.x[f] for f in model.f) >= model.b[n]
model.nb = Constraint(model.n, rule=nutrient_balance_rule, doc='nutrient balance (units)')

def cost_balance_rule(model):
    return model.cost == sum(model.x[f] for f in model.f)
model.cb = Constraint(rule=cost_balance_rule, doc='cost balance (dollars)')
model.obj = Objective(expr=model.cost, sense=minimize, doc='objective function')
PYEOF
```

</details>

<details>
<summary>diet_data.json (click to expand)</summary>

```bash
cat > /tmp/optiskill_workspace/diet_data.json << 'JSONEOF'
{
  "nutrient_requirements": {
    "calorie": 3, "protein": 70, "calcium": 0.8, "iron": 12,
    "vitamin-a": 5, "vitamin-b1": 1.8, "vitamin-b2": 2.7,
    "niacin": 18, "vitamin-c": 75
  },
  "nutritive_values": {
    "wheat":     {"calorie":44.7,"protein":1411,"calcium":2.0,"iron":365,"vitamin-a":0,"vitamin-b1":55.4,"vitamin-b2":33.3,"niacin":441,"vitamin-c":0},
    "cornmeal":  {"calorie":36,"protein":897,"calcium":1.7,"iron":99,"vitamin-a":30.9,"vitamin-b1":17.4,"vitamin-b2":7.9,"niacin":106,"vitamin-c":0},
    "cannedmilk":{"calorie":8.4,"protein":422,"calcium":15.1,"iron":9,"vitamin-a":26,"vitamin-b1":3,"vitamin-b2":23.5,"niacin":11,"vitamin-c":60},
    "margarine": {"calorie":20.6,"protein":17,"calcium":0.6,"iron":6,"vitamin-a":55.8,"vitamin-b1":0.2,"vitamin-b2":0,"niacin":0,"vitamin-c":0},
    "cheese":    {"calorie":7.4,"protein":448,"calcium":16.4,"iron":19,"vitamin-a":28.1,"vitamin-b1":0.8,"vitamin-b2":10.3,"niacin":4,"vitamin-c":0},
    "peanut-b":  {"calorie":15.7,"protein":661,"calcium":1,"iron":48,"vitamin-a":0,"vitamin-b1":9.6,"vitamin-b2":8.1,"niacin":471,"vitamin-c":0},
    "lard":      {"calorie":41.7,"protein":0,"calcium":0,"iron":0,"vitamin-a":0,"vitamin-b1":0,"vitamin-b2":0.5,"niacin":5,"vitamin-c":0},
    "liver":     {"calorie":2.2,"protein":333,"calcium":0.2,"iron":139,"vitamin-a":169.2,"vitamin-b1":6.4,"vitamin-b2":50.8,"niacin":316,"vitamin-c":525},
    "porkroast": {"calorie":4.4,"protein":249,"calcium":0.3,"iron":37,"vitamin-a":0,"vitamin-b1":18.2,"vitamin-b2":3.6,"niacin":79,"vitamin-c":0},
    "salmon":    {"calorie":5.8,"protein":705,"calcium":6.8,"iron":45,"vitamin-a":3.5,"vitamin-b1":1,"vitamin-b2":4.9,"niacin":209,"vitamin-c":0},
    "greenbeans":{"calorie":2.4,"protein":138,"calcium":3.7,"iron":80,"vitamin-a":69,"vitamin-b1":4.3,"vitamin-b2":5.8,"niacin":37,"vitamin-c":862},
    "cabbage":   {"calorie":2.6,"protein":125,"calcium":4,"iron":36,"vitamin-a":7.2,"vitamin-b1":9,"vitamin-b2":4.5,"niacin":26,"vitamin-c":5369},
    "onions":    {"calorie":5.8,"protein":166,"calcium":3.8,"iron":59,"vitamin-a":16.6,"vitamin-b1":4.7,"vitamin-b2":5.9,"niacin":21,"vitamin-c":1184},
    "potatoes":  {"calorie":14.3,"protein":336,"calcium":1.8,"iron":118,"vitamin-a":6.7,"vitamin-b1":29.4,"vitamin-b2":7.1,"niacin":198,"vitamin-c":2522},
    "spinach":   {"calorie":1.1,"protein":106,"calcium":0,"iron":138,"vitamin-a":918.4,"vitamin-b1":5.7,"vitamin-b2":13.8,"niacin":33,"vitamin-c":2755},
    "sweet-pot": {"calorie":9.6,"protein":138,"calcium":2.7,"iron":54,"vitamin-a":290.7,"vitamin-b1":8.4,"vitamin-b2":5.4,"niacin":83,"vitamin-c":1912},
    "peaches":   {"calorie":8.5,"protein":87,"calcium":1.7,"iron":173,"vitamin-a":86.8,"vitamin-b1":1.2,"vitamin-b2":4.3,"niacin":55,"vitamin-c":57},
    "prunes":    {"calorie":12.8,"protein":99,"calcium":2.5,"iron":154,"vitamin-a":85.7,"vitamin-b1":3.9,"vitamin-b2":4.3,"niacin":65,"vitamin-c":257},
    "limabeans": {"calorie":17.4,"protein":1055,"calcium":3.7,"iron":459,"vitamin-a":5.1,"vitamin-b1":26.9,"vitamin-b2":38.2,"niacin":93,"vitamin-c":0},
    "navybeans": {"calorie":26.9,"protein":1691,"calcium":11.4,"iron":792,"vitamin-a":0,"vitamin-b1":38.4,"vitamin-b2":24.6,"niacin":217,"vitamin-c":0}
  }
}
JSONEOF
```

</details>

### 0c. Build, solve, and store the base model

```bash
python3 << 'BOOTSTRAP'
import json, os, sys, copy
import pyomo.environ as pe
from pyomo.environ import value, SolverFactory, Var, Param, Set, Constraint, Objective, Suffix
import cloudpickle, pickle

WORKSPACE = "/tmp/optiskill_workspace"

# ── Helper: solver fallback chain ──
def get_solver():
    """Return the best available solver (Gurobi > HiGHS > GLPK)."""
    from pyomo.opt import SolverFactory
    for name in ["gurobi", "appsi_highs", "glpk"]:
        try:
            s = SolverFactory(name)
            if s.available():
                print(f"Using solver: {name}")
                return s
        except Exception:
            continue
    raise RuntimeError("No solver found. Install highspy or glpk.")

# ── Load data ──
with open(f"{WORKSPACE}/diet_data.json") as f:
    data = json.load(f)

# ── Execute model source ──
namespace = {"data": data}
with open(f"{WORKSPACE}/diet.py") as f:
    exec(f.read(), namespace)
model = namespace["model"]

# ── Add dual suffix for sensitivity analysis ──
model.dual = Suffix(direction=Suffix.IMPORT)

# ── Solve ──
solver = get_solver()
results = solver.solve(model, tee=False)
status = str(results.solver.termination_condition)
print(f"Solver status: {status}")
print(f"Objective:     {value(model.obj):.6f}")

# ── Save pickle ──
pkl_path = f"{WORKSPACE}/diet.pkl"
with open(pkl_path, "wb") as f:
    cloudpickle.dump(model, f, protocol=pickle.HIGHEST_PROTOCOL)

# ── Extract model info ──
def extract_model_info(m, solver_status):
    info = {"params": {}, "vars": {}, "constraints": {}, "sets": {}, "obj": {}}
    for obj in m.component_objects(Objective, active=True):
        info["obj"] = {"name": obj.name, "sense": str(obj.sense),
                       "value": value(obj), "sol_status": solver_status}
    for s in m.component_objects(Set, active=True):
        info["sets"][s.name] = {"doc": s.doc or "", "members": list(s)}
    for p in m.component_objects(Param, active=True):
        pd = {}
        for idx in p:
            k = str(idx) if idx is not None else "_scalar"
            pd[k] = {"value": value(p[idx]), "mutable": p._mutable}
        info["params"][p.name] = {"doc": p.doc or "", "data": pd}
    for v in m.component_objects(Var, active=True):
        vd = {}
        for idx in v:
            k = str(idx) if idx is not None else "_scalar"
            try: sol = value(v[idx])
            except: sol = None
            vd[k] = {"solution": sol, "lb": v[idx].lb, "ub": v[idx].ub}
        info["vars"][v.name] = {"doc": v.doc or "", "data": vd}
    for c in m.component_objects(Constraint, active=True):
        cd = {}
        for idx in c:
            k = str(idx) if idx is not None else "_scalar"
            try:
                bv = value(c[idx].body)
                lb = float(c[idx].lower) if c[idx].lower is not None else None
                ub = float(c[idx].upper) if c[idx].upper is not None else None
                slack = min((bv - lb) if lb is not None else 1e99,
                            (ub - bv) if ub is not None else 1e99)
                binding = abs(slack) < 1e-6
            except:
                bv, lb, ub, binding = None, None, None, None
            cd[k] = {"body_value": bv, "lower": lb, "upper": ub, "is_binding": binding}
        info["constraints"][c.name] = {"doc": c.doc or "", "data": cd}
    return info

model_info = extract_model_info(model, status)

# ── Store JSON snapshot ──
import json as _json
with open(f"{WORKSPACE}/models_dict.json", "w") as f:
    safe = {}
    for kk, vv in {"diet": {**model_info, "local_path_to_object": pkl_path}}.items():
        safe[kk] = {k2: v2 for k2, v2 in vv.items() if k2 != "model_object"}
    _json.dump(safe, f, default=str, indent=2)

print(f"\nBootstrap complete. Models: ['diet']")
print(f"Params: {list(model_info['params'].keys())}")
print(f"Vars:   {list(model_info['vars'].keys())}")
print(f"Constrs:{list(model_info['constraints'].keys())}")
BOOTSTRAP
```

**Expected output:**
```
Solver status: optimal
Objective:     0.108662
Bootstrap complete. Models: ['diet']
Params: ['b', 'a']
Vars:   ['x', 'cost']
Constrs:['nb', 'cb']
```

---

## Query Classification

Classify each user question and jump to the matching workflow section:

```
1. "What does this model do?" / "Explain the model"     → Workflow A: Model Description
2. "What is the value of X?" / "Which constraints bind?" → Workflow B: Retrieval
3. "How sensitive is cost to X?" / "Shadow prices?"      → Workflow C: Sensitivity
4. "What if X increases by 20%?" / "Set X to 500"        → Workflow D: What-If
5. "Why doesn't the model choose X?" / "Force X on"      → Workflow E: Why-Not
6. "Add constraint..." / "Compute..." / anything else     → Workflow F: Code Execution
```

**Key distinction for D vs E:**
- Target is a **Param** (input data) → **What-If** (Workflow D)
- Target is a **Var** (decision) → **Why-Not** (Workflow E)

**Key distinction for C vs D:**
- No specific value given (direction/rate) → **Sensitivity** (Workflow C)
- Specific value given ("by 20%", "set to 500") → **What-If** (Workflow D)

---

## Shared Utilities

Every workflow below uses these helpers. They are defined in the bootstrap and
should be re-declared in each heredoc for self-containedness:

```python
# ── Solver fallback ──
def get_solver():
    from pyomo.opt import SolverFactory
    for name in ["gurobi", "appsi_highs", "glpk"]:
        try:
            s = SolverFactory(name)
            if s.available():
                return s
        except Exception:
            continue
    raise RuntimeError("No solver found. Install highspy or glpk.")

# ── CRITICAL: always use value() to read Pyomo objects ──
from pyomo.environ import value  # value(model.x[i]) not float(model.x[i])
```

---

## Workflow A — Model Description

**Goal:** Generate a narrative overview and structured component reference.

### A1. Load and summarise components

```bash
python3 << 'WF_A'
import json, cloudpickle
from pyomo.environ import value

WORKSPACE = "/tmp/optiskill_workspace"

with open(f"{WORKSPACE}/models_dict.json") as f:
    models_dict = json.load(f)
with open(f"{WORKSPACE}/diet.pkl", "rb") as f:
    model = cloudpickle.load(f)

info = models_dict["diet"]

# ── Part 1: Narrative ──
print("=" * 60)
print("PART 1: MODEL NARRATIVE")
print("=" * 60)
print(f"\nObjective: {info['obj']['name']} ({info['obj']['sense']})")
print(f"Status: {info['obj']['sol_status']}, Value: {info['obj']['value']:.6f}\n")

print("SETS:")
for sname, sinfo in info["sets"].items():
    members = sinfo["members"]
    sample = ", ".join(str(m) for m in members[:5])
    print(f"  {sname}: {sinfo['doc']} — {len(members)} items ({sample}{'...' if len(members) > 5 else ''})")

print("\nPARAMETERS:")
for pname, pinfo in info["params"].items():
    print(f"  {pname}: {pinfo['doc']} ({len(pinfo['data'])} entries, mutable={list(pinfo['data'].values())[0]['mutable']})")

print("\nVARIABLES:")
for vname, vinfo in info["vars"].items():
    print(f"  {vname}: {vinfo['doc']}")
    nonzero = [(idx, v["solution"]) for idx, v in vinfo["data"].items()
               if v["solution"] is not None and abs(v["solution"]) > 1e-6]
    for idx, sol in nonzero[:3]:
        print(f"    [{idx}] = {sol:.6f}")

print("\nCONSTRAINTS:")
for cname, cinfo in info["constraints"].items():
    binding = [idx for idx, cd in cinfo["data"].items() if cd.get("is_binding")]
    print(f"  {cname}: {cinfo['doc']} — {len(binding)}/{len(cinfo['data'])} binding")

# ── Part 2: Component Reference ──
print("\n" + "=" * 60)
print("PART 2: COMPONENT REFERENCE")
print("=" * 60)

print("\nPARAMETER ROLES:")
for pname, pinfo in info["params"].items():
    print(f"  {pname} | {pinfo['doc']}")

print("\nVARIABLE DOMAINS:")
for vname, vinfo in info["vars"].items():
    sample = list(vinfo["data"].values())[0]
    print(f"  {vname} | {vinfo['doc']} | bounds=[{sample['lb']}, {sample['ub']}]")

print("\nCONSTRAINT STRUCTURE:")
for cname, cinfo in info["constraints"].items():
    sample = list(cinfo["data"].values())[0]
    print(f"  {cname} | {cinfo['doc']} | lower={sample.get('lower')}, upper={sample.get('upper')}")

print("\nUse Part 2 names exactly when referencing components in other workflows.")
WF_A
```

**Expected output:** A two-part description with narrative overview and component reference table listing all sets, parameters, variables, constraints, and objective with their Pyomo names.

---

## Workflow B — Retrieval

**Goal:** Read current values from the solved model without modification.

### B1. Query specific components

```bash
python3 << 'WF_B'
import json, cloudpickle
from pyomo.environ import value

WORKSPACE = "/tmp/optiskill_workspace"

with open(f"{WORKSPACE}/diet.pkl", "rb") as f:
    model = cloudpickle.load(f)
with open(f"{WORKSPACE}/models_dict.json") as f:
    info = json.load(f)["diet"]

# ── Retrieve parameter values ──
print("=== Parameter b (nutrient requirements) ===")
for nutrient in model.b:
    print(f"  b['{nutrient}'] = {value(model.b[nutrient])}")

# ── Retrieve variable values (non-zero only) ──
print("\n=== Variable x (food spending, non-zero) ===")
for food in model.f:
    amt = value(model.x[food])
    if amt > 1e-6:
        print(f"  x['{food}'] = ${amt:.6f}/day")

# ── Retrieve constraint status ──
print("\n=== Constraint Status ===")
for cname, cinfo in info["constraints"].items():
    print(f"\n{cname}: {cinfo['doc']}")
    for idx, cd in list(cinfo["data"].items())[:5]:
        binding = "[BINDING]" if cd.get("is_binding") else "[slack]"
        print(f"  [{idx}] {binding} body={cd['body_value']:.4f}")

# ── Objective ──
print(f"\n=== Objective ===")
print(f"  {info['obj']['name']}: {info['obj']['value']:.6f} ({info['obj']['sol_status']})")
WF_B
```

**Expected output:** Lists of parameter values, non-zero variable solutions, constraint binding status, and objective value. Adapt the component names and indices to the user's specific query.

---

## Workflow C — Sensitivity Analysis (Dual Values)

**Goal:** Quantify marginal impact of parameter changes using shadow prices.

### C1. Check for integer variables (duals require LP)

```bash
python3 << 'WF_C1'
import cloudpickle
from pyomo.environ import Var

WORKSPACE = "/tmp/optiskill_workspace"
with open(f"{WORKSPACE}/diet.pkl", "rb") as f:
    model = cloudpickle.load(f)

has_integer = any(v[idx].is_integer() or v[idx].is_binary()
                  for v in model.component_objects(Var, active=True)
                  for idx in v)

if has_integer:
    print("WARNING: MILP detected — dual values do not apply.")
    print("Use Workflow D (What-If) with specific parameter values instead.")
else:
    print("LP model confirmed. Dual values are applicable.")
WF_C1
```

### C2. Extract and interpret dual values

```bash
python3 << 'WF_C2'
import cloudpickle, copy
from pyomo.environ import value, Suffix, SolverFactory

WORKSPACE = "/tmp/optiskill_workspace"
with open(f"{WORKSPACE}/diet.pkl", "rb") as f:
    base_model = cloudpickle.load(f)

model = copy.deepcopy(base_model)

# Ensure dual suffix exists
if not hasattr(model, 'dual'):
    model.dual = Suffix(direction=Suffix.IMPORT)

def get_solver():
    for name in ["gurobi", "appsi_highs", "glpk"]:
        try:
            s = SolverFactory(name)
            if s.available():
                return s
        except Exception:
            continue
    raise RuntimeError("No solver found.")

solver = get_solver()
results = solver.solve(model, tee=False)
print(f"Solver: {results.solver.termination_condition}\n")

print("=== SENSITIVITY ANALYSIS: Dual Values (Shadow Prices) ===\n")

# Read duals for nutrient balance constraints
ranking = []
for nutrient in model.n:
    dual_val = model.dual.get(model.nb[nutrient], 0)
    current_req = value(model.b[nutrient])
    intake = sum(value(model.a[f, nutrient]) * value(model.x[f]) for f in model.f)
    slack = intake - current_req
    binding = "BINDING" if abs(slack) < 1e-6 else "slack"

    print(f"{nutrient:15s}: dual={dual_val:10.6f}  ({binding})")
    if abs(dual_val) > 1e-6:
        print(f"  -> +1 unit requirement costs ${abs(dual_val):.6f} more")
    ranking.append((nutrient, abs(dual_val)))

print("\n=== SENSITIVITY RANKING (highest impact first) ===")
for nutrient, d in sorted(ranking, key=lambda x: -x[1])[:5]:
    print(f"  {nutrient}: ${d:.6f}/unit")
WF_C2
```

**Expected output:** For each constraint, the dual value, binding status, and marginal cost interpretation. A ranking identifies the highest-leverage parameters.

---

## Workflow D — What-If Analysis

**Goal:** Modify a parameter to a specific value, re-solve, and compare.

### D1. Modify parameter and re-solve

```bash
python3 << 'WF_D'
import json, copy, cloudpickle, pickle
from pyomo.environ import value, SolverFactory

WORKSPACE = "/tmp/optiskill_workspace"

with open(f"{WORKSPACE}/diet.pkl", "rb") as f:
    base_model = cloudpickle.load(f)

base_obj = value(base_model.obj)
model = copy.deepcopy(base_model)

# === USER SCENARIO: Double the protein requirement ===
# Mapping: "increase by X%" → param[idx] = value(param[idx]) * (1 + X/100)
#          "set to X"        → param[idx] = X
#          "make unavailable" → param[idx] = 0
TARGET_PARAM = "b"
TARGET_INDEX = "protein"
old_val = value(model.b[TARGET_INDEX])
new_val = old_val * 2  # Example: double
model.b[TARGET_INDEX] = new_val
print(f"Changed {TARGET_PARAM}['{TARGET_INDEX}']: {old_val} -> {new_val}")

# Re-solve
def get_solver():
    for name in ["gurobi", "appsi_highs", "glpk"]:
        try:
            s = SolverFactory(name)
            if s.available():
                return s
        except Exception:
            continue
    raise RuntimeError("No solver found.")

solver = get_solver()
results = solver.solve(model, tee=False)
status = str(results.solver.termination_condition)
new_obj = value(model.obj)

print(f"\nStatus: {status}")
if status == "optimal":
    delta = new_obj - base_obj
    pct = 100 * delta / base_obj if base_obj else float('nan')
    print(f"Objective: {base_obj:.6f} -> {new_obj:.6f} ({delta:+.6f}, {pct:+.1f}%)")

    print("\nVariable changes:")
    for f_idx in model.f:
        old_x = value(base_model.x[f_idx])
        new_x = value(model.x[f_idx])
        if abs(new_x - old_x) > 1e-6:
            print(f"  x['{f_idx}']: {old_x:.4f} -> {new_x:.4f}")

    # Save modified model
    version = f"diet__{TARGET_PARAM}_{TARGET_INDEX}_modified"
    with open(f"{WORKSPACE}/{version}.pkl", "wb") as f:
        cloudpickle.dump(model, f, protocol=pickle.HIGHEST_PROTOCOL)
    print(f"\nSaved as '{version}'")
else:
    print(f"Scenario made the model {status}.")
    print("Investigate constraint conflicts or try a less extreme change.")
WF_D
```

**Expected output:** Old vs. new objective, percentage change, and which variables shifted. Adapt `TARGET_PARAM`, `TARGET_INDEX`, and the modification formula to the user's query.

---

## Workflow E — Why-Not Analysis

**Goal:** Force a decision variable to a user-specified value, re-solve, and explain the cost penalty.

### E1. Force variable and compare

```bash
python3 << 'WF_E'
import json, copy, cloudpickle, pickle
from pyomo.environ import value, SolverFactory

WORKSPACE = "/tmp/optiskill_workspace"

with open(f"{WORKSPACE}/diet.pkl", "rb") as f:
    base_model = cloudpickle.load(f)

base_obj = value(base_model.obj)
model = copy.deepcopy(base_model)

# === USER SCENARIO: Force liver spending to $0 (ban liver) ===
TARGET_VAR = "x"
TARGET_INDEX = "liver"
FORCED_VALUE = 0.0

old_val = value(model.x[TARGET_INDEX])
model.x[TARGET_INDEX].fix(FORCED_VALUE)
print(f"Forced {TARGET_VAR}['{TARGET_INDEX}'] = {FORCED_VALUE} (was {old_val:.6f})")

# Re-solve
def get_solver():
    for name in ["gurobi", "appsi_highs", "glpk"]:
        try:
            s = SolverFactory(name)
            if s.available():
                return s
        except Exception:
            continue
    raise RuntimeError("No solver found.")

solver = get_solver()
results = solver.solve(model, tee=False)
status = str(results.solver.termination_condition)
new_obj = value(model.obj)

print(f"\nStatus: {status}")
if status == "optimal":
    delta = new_obj - base_obj
    pct = 100 * delta / base_obj if base_obj else float('nan')
    print(f"Objective: {base_obj:.6f} -> {new_obj:.6f} ({delta:+.6f}, {pct:+.1f}%)")

    if abs(delta) < 1e-8:
        print("\nNo cost penalty — this decision was already optimal or irrelevant.")
    else:
        print(f"\nCost penalty of forcing this decision: ${abs(delta):.6f}/day")

    print("\nSubstitution effects:")
    for f_idx in model.f:
        old_x = value(base_model.x[f_idx])
        new_x = value(model.x[f_idx])
        if abs(new_x - old_x) > 1e-6:
            print(f"  x['{f_idx}']: {old_x:.4f} -> {new_x:.4f} ({new_x - old_x:+.4f})")

    # Save forced model
    version = f"diet__{TARGET_VAR}_{TARGET_INDEX}_forced"
    with open(f"{WORKSPACE}/{version}.pkl", "wb") as f:
        cloudpickle.dump(model, f, protocol=pickle.HIGHEST_PROTOCOL)
    print(f"\nSaved as '{version}'")
else:
    print(f"Forcing this decision made the model {status}.")
    print("The forced value conflicts with existing constraints.")
WF_E
```

**Expected output:** Cost penalty, substitution effects showing which other variables shifted. Adapt `TARGET_VAR`, `TARGET_INDEX`, and `FORCED_VALUE` to the user's query. For binary variables, use `.fix(1)` (force on) or `.fix(0)` (force off).

---

## Workflow F — Python Code Execution

**Goal:** Execute custom Python code for complex model manipulations.

### F1. Common patterns

```bash
python3 << 'WF_F'
import json, copy, cloudpickle
from pyomo.environ import value, SolverFactory, Constraint, Suffix

WORKSPACE = "/tmp/optiskill_workspace"

with open(f"{WORKSPACE}/diet.pkl", "rb") as f:
    base_model = cloudpickle.load(f)

model = copy.deepcopy(base_model)
base_obj = value(model.obj)

# ── PATTERN A: Add a new constraint ──
min_liver = 0.01
model.liver_min = Constraint(expr=model.x['liver'] >= min_liver)
print(f"Added constraint: x['liver'] >= {min_liver}")

# ── PATTERN B: Deactivate a constraint ──
# model.nb['calorie'].deactivate()
# print("Deactivated calorie constraint")

# ── Re-solve ──
def get_solver():
    for name in ["gurobi", "appsi_highs", "glpk"]:
        try:
            s = SolverFactory(name)
            if s.available():
                return s
        except Exception:
            continue
    raise RuntimeError("No solver found.")

solver = get_solver()
results = solver.solve(model, tee=False)
status = str(results.solver.termination_condition)
new_obj = value(model.obj)

print(f"\nStatus: {status}")
print(f"Objective: {base_obj:.6f} -> {new_obj:.6f} ({new_obj - base_obj:+.6f})")

# ── PATTERN C: Extract dual values ──
model2 = copy.deepcopy(base_model)
model2.dual = Suffix(direction=Suffix.IMPORT)
solver.solve(model2, tee=False)
print("\nDual values (nutrient constraints):")
for n in model2.n:
    d = model2.dual.get(model2.nb[n], 0)
    print(f"  nb['{n}']: {d:.6f}")

# ── PATTERN D: Iterative cost breakdown ──
print("\nCost breakdown:")
total = sum(value(base_model.x[f]) for f in base_model.f)
for f in base_model.f:
    amt = value(base_model.x[f])
    if amt > 1e-6:
        print(f"  {f}: ${amt:.6f} ({100*amt/total:.1f}%)")
WF_F
```

**Expected output:** Results of the custom code execution. Adapt the patterns to the user's specific request. Always use `value()` for Pyomo objects and `copy.deepcopy()` to preserve the base model.

### Critical Pyomo Rules

```python
# CORRECT                                    # WRONG
old = value(model.b['protein'])              # old = float(model.b['protein'])  # TypeError
model.b['protein'] = old * 2                 # model.b['protein'] *= 2          # TypeError
total = sum(value(model.x[f]) for f in model.f)  # sum(model.x[f] ...)         # Symbolic, not numeric
```

---

## Generalisability

This skill applies to **any** Pyomo ConcreteModel, not just the bundled diet example.

To adapt to a different model:

1. Replace `diet.py` and `diet_data.json` in Prerequisites 0b with your model source and data.
2. All six workflows use the same pattern: load pickle → deep copy → modify → re-solve → compare.
3. Component names (`model.b`, `model.x`, `model.nb`) change per model — use Workflow A to discover them.
4. The same approach generalises to other algebraic modelling languages (AMPL, GAMS, JuMP) by substituting the equivalent component-access API.

**Solver compatibility:** Workflows use Gurobi if licensed, otherwise HiGHS (LP/MILP/QP) or GLPK as fallback. No commercial licence required for the open-source path.

---

## Edge Cases

- **MILP models**: Dual values (Workflow C) do not apply. Use Workflow D (What-If) for concrete scenario comparison instead.
- **Infeasible scenarios**: If a What-If or Why-Not re-solve returns "infeasible", report this clearly and suggest the user try a less extreme parameter change.
- **Non-mutable parameters**: If `mutable=False`, Pyomo blocks runtime changes. Inform the user that the parameter is fixed by design.
- **Large index sets**: Filter output to show only significant values (e.g., non-zero variables, binding constraints).
- **Unsolved models**: If `sol_status == "unknown"`, solve the model first before retrieving values.
- **Multiple simultaneous changes**: In What-If, modify all parameters before the single `solver.solve()` call.
- **Binary variables in Why-Not**: Use `.fix(1)` for "must open" or `.fix(0)` for "must close".