Updated matrix F in Palma Ratio Linearization

CLAUDE.md

@@ -0,0 +1,90 @@
+# CLAUDE.md
+
+This file provides guidance to Claude Code (claude.ai/code) when working with code in this repository.
+
+## Project Overview
+
+FairLoadDelivery is a Julia package for fair load shedding optimization in power distribution systems. It implements a bilevel optimization framework that balances efficiency with fairness metrics when determining which loads to shed during capacity shortages.
+
+## Build & Development Commands
+
+```bash
+# Activate and install dependencies
+julia --project=. -e 'using Pkg; Pkg.instantiate()'
+
+# Run tests
+julia --project=. -e 'using Pkg; Pkg.test()'
+
+# Start interactive development session
+julia --project=.
+
+# In Julia REPL:
+using Revise
+using FairLoadDelivery
+```
+
+## Running Scripts
+
+```julia
+# Main FLDP workflow
+include("script/FLDP.jl")
+
+# Network setup example
+eng, math, lbs, critical_id = setup_network("ieee_13_aw_edit/motivation_b.dss", 0.9, [])
+
+# Solve MLD problems
+pm_soln = FairLoadDelivery.solve_mc_mld_switch_integer(math, gurobi)
+pm_soln = solve_mc_mld_shed_implicit_diff(math, ipopt)
+```
+
+## Architecture
+
+### Bilevel Optimization Framework
+
+- **Upper Level**: Optimizes fairness metrics by adjusting load weights
+- **Lower Level**: Solves Minimum Load Delivery (MLD) problem given weights, using implicit differentiation (`DiffOpt`) to compute gradients through the optimization
+
+### Core Components
+
+- `src/core/` - Variable definitions, constraints, and objectives for JuMP models
+- `src/prob/` - Problem formulations (OPF, MLD, Power Flow)
+- `src/implementation/` - Algorithm implementations:
+  - `network_setup.jl` - Parses OpenDSS files, identifies load blocks
+  - `lower_level_mld.jl` - Lower-level solver with DiffOpt differentiation
+  - `palma_relaxation.jl` - Palma ratio fairness optimization
+  - `random_rounding.jl` - Converts relaxed solutions to integer-feasible
+  - `other_fair_funcs.jl` - Fairness metrics (Jain, proportional, min-max, efficiency)
+
+### Key Concepts
+
+- **Load Blocks**: Connected regions of loads that can be shed together, enforced by switch topology
+- **Three-Phase Unbalanced**: Uses PowerModelsDistribution for realistic distribution system modeling
+- **Radiality Constraints**: Ensures tree topology in distribution networks
+
+### Solvers
+
+Preconfigured optimizers are exported from the module:
+- `ipopt` - NLP/continuous problems (primary)
+- `gurobi` - Mixed-integer programming
+- `highs` - Secondary MIP solver
+- `juniper` - Mixed-integer nonlinear
+
+### Data Format
+
+Network files use OpenDSS `.dss` format, located in `data/`. Primary test case: `ieee_13_aw_edit/motivation_b.dss`
+
+## Key Patterns
+
+### Reference Extensions
+
+Network preprocessing uses the extension pattern:
+```julia
+ref_extensions=[ref_add_load_blocks!]
+# or for rounded solutions:
+ref_extensions=[ref_add_rounded_load_blocks!]
+```
+
+### Constants
+
+- `zero_tol = 1e-9` - Numerical tolerance for zero checks
+- Module aliases: `_PMD` for PowerModelsDistribution, `_IM` for InfrastructureModels

script/reformulation/README.md

@@ -0,0 +1,176 @@
+# Palma Ratio Reformulation
+
+This directory contains the reformulated Palma ratio optimization for fair load shedding.
+
+## Overview
+
+The Palma ratio measures inequality as:
+```
+Palma = sum(top 10% load shed) / sum(bottom 40% load shed)
+```
+
+Lower Palma ratio = more equal distribution of load shedding.
+
+## Files
+
+| File | Description |
+|------|-------------|
+| `load_shed_as_parameter.jl` | Core optimization: Palma ratio minimization with McCormick envelopes and Charnes-Cooper transformation |
+| `opendss_experiment.jl` | Bilevel optimization loop integrating with OpenDSS network data via DiffOpt |
+| `validate_reformulation.jl` | Unit tests for the reformulation |
+| `diagnose_jacobian.jl` | Diagnostic script to verify Jacobian computation |
+
+## Mathematical Formulation
+
+### Decision Variables
+- `Δw[j]`: Weight changes (continuous, bounded by trust region)
+- `a[i,j]`: Permutation matrix (binary or relaxed)
+- `u[i,j]`: McCormick auxiliary for `a[i,j] * pshed_new[j]`
+- `σ`: Charnes-Cooper scaling variable
+
+### Key Expressions
+```julia
+# P_shed as Taylor expansion (NOT a variable)
+pshed_new[j] = pshed_prev[j] + Σₖ (∂pshed[j]/∂w[k]) * Δw[k]
+
+# Sorted values via permutation
+sorted[i] = Σⱼ u[i,j]  where u[i,j] ≈ a[i,j] * pshed_new[j]
+```
+
+### Charnes-Cooper Transformation
+Converts ratio minimization to linear objective:
+```
+min (top 10% sum) / (bottom 40% sum)
+→  min (top 10% sum) * σ   s.t. (bottom 40% sum) * σ = 1
+```
+
+## Known Issues and Limitations
+
+### 1. McCormick Relaxation vs Binary Permutation
+
+**Problem**: McCormick envelopes are tight only at binary (0,1) points.
+
+| Setting | Pros | Cons |
+|---------|------|------|
+| `relax_binary=true` | Fast LP solve | **Sorting breaks** - u values collapse to zero |
+| `relax_binary=false` | Correct sorting | **Slow MIQP** - 225 binary variables for n=15 |
+
+**Current approach**: Use `relax_binary=false` with Gurobi time limits.
+
+### 2. Taylor Approximation Validity
+
+The bilevel optimization relies on:
+```
+pshed_new ≈ pshed_prev + Jacobian * Δw
+```
+
+This is only valid for **small** weight changes. If:
+- Trust region is too large
+- System is highly nonlinear
+- Jacobian is inaccurate
+
+The predicted Palma ratio may be **very different** from actual.
+
+### 3. Jacobian Computation (CRITICAL BUG)
+
+**Status: DiffOpt sensitivities are BROKEN for this NLP**
+
+The Jacobian `∂pshed/∂w` computed via DiffOpt forward differentiation produces **WRONG** values:
+
+| Issue | Expected | Actual (DiffOpt) |
+|-------|----------|------------------|
+| Sign | Negative | **Positive** (all 15 loads) |
+| Magnitude | O(demand) ≈ 17-1155 | **O(100,000)** - 100x too large |
+
+**Root cause**: DiffOpt's KKT-based differentiation fails on this non-convex NLP because:
+1. "Inertia correction needed" warnings indicate ill-conditioned KKT matrix
+2. Non-convex power flow constraints create saddle points
+3. Relaxed binary indicators (z_demand) cause degeneracy
+
+**Workaround**: Use **finite differences** instead of DiffOpt:
+```julia
+# Actually perturb weights and re-solve MLD
+δ = 0.01
+for j in 1:n_loads
+    math_perturbed = deepcopy(math)
+    math_perturbed["load"][j]["weight"] += δ
+    # ... re-solve and compute (pshed_new - pshed_baseline) / δ
+end
+```
+
+The `diagnose_jacobian.jl` script now includes both DiffOpt and finite difference comparison.
+Run it to verify:
+```bash
+julia --project=. script/reformulation/diagnose_jacobian.jl
+```
+
+### 4. Charnes-Cooper Creates Quadratic Constraints
+
+The constraint `(bottom 40% sum) * σ = 1` is bilinear, requiring:
+- Gurobi with `NonConvex=2`, or
+- Ipopt (NLP solver)
+
+HiGHS cannot be used (LP/MILP only).
+
+## Alternative Approaches
+
+Located in `src/implementation/other_fair_funcs.jl`:
+
+| Metric | Formula | Complexity | Recommendation |
+|--------|---------|------------|----------------|
+| **Proportional Fairness** | `max Σ log(pshed + ε)` | NLP | ⭐ Best alternative |
+| **Jain's Index** | `max (Σpshed)² / (n·Σpshed²)` | NLP | Good for bounded metric |
+| **Min-Max** | `max min(pshed)` | LP | Extreme fairness |
+| **Palma Ratio** | `min top10% / bottom40%` | MIQP | Current (slow) |
+
+**Recommendation**: If Palma MIQP is too slow, use **Proportional Fairness** - it has similar fairness properties without requiring sorting/permutation matrices.
+
+## Usage
+
+### Basic experiment
+```julia
+include("script/reformulation/opendss_experiment.jl")
+run_mwe(ls_percent=0.5)
+```
+
+### With custom settings
+```julia
+result = solve_palma_ratio_minimization(
+    "ieee_13_aw_edit/motivation_a.dss";
+    ls_percent = 0.5,        # 50% capacity → forces load shedding
+    iterations = 5,
+    trust_radius = 0.1,      # Smaller = more conservative
+    experiment_name = "my_experiment"
+)
+```
+
+### Diagnose Jacobian issues
+```bash
+julia --project=. script/reformulation/diagnose_jacobian.jl
+```
+
+## Debugging
+
+### Palma ratio exploding?
+1. Check Jacobian diagonal: should be **negative**
+2. Check predicted vs actual Palma after each iteration
+3. Reduce trust radius (try 0.01 instead of 0.1)
+4. Check if bottom 40% sum → 0 (makes Palma undefined)
+
+### Gurobi timing out?
+1. Increase `MIPGap` to 0.01 (1%)
+2. Reduce time limit
+3. Consider switching to Proportional Fairness
+
+### Segmentation fault?
+Usually caused by Julia/solver memory issues. Try:
+- Restart Julia session
+- Reduce number of iterations
+- Use fewer loads
+
+## References
+
+- Charnes-Cooper transformation: Charnes & Cooper (1962)
+- McCormick envelopes: McCormick (1976)
+- Palma ratio: Cobham & Sumner (2013)
+- DiffOpt.jl: https://github.com/jump-dev/DiffOpt.jl