Skip to content

Conversation

@samtalki
Copy link
Collaborator

@samtalki samtalki commented Jan 5, 2026

Matrix F Corrections for Palma Ratio Linearization

This PR suggests some changes to the constraint matrix F and related components in the Palma ratio linearization implementation. It attempts to resolve the bug mentioned in #3

I tried to align things with the mathematical formulation on the Overleaf. I used Claude Code Pro with the Opus 4.5 model to try to help.

Possible issues

I believe there are three possible issues in the current implementation:

  1. Indexing Inconsistency: Mixed row-major and column-major vectorization conventions
  2. McCormick Envelope RHS Errors: Possibly incorrect right-hand side values for McCormick constraints
  3. Undefined Variable References: A_long used undefined variable A

Issue 1: Indexing Inconsistency

Problem

The original code mixed two different vectorization conventions for the permutation matrix:

Component Original Indexing
A, AT, T Row-major: k = (i-1)*n + j
A_row, A_col Column-major: k = i + (j-1)*n

This caused the sorting constraint (using T) to operate inconsistently with the permutation constraints (using A_row, A_col).

Example for n=3

For a permutation matrix entry a_{i,j}:

  • Row-major: vec(a) = [a_{11}, a_{12}, a_{13}, a_{21}, a_{22}, a_{23}, a_{31}, a_{32}, a_{33}]
  • Column-major: vec(a) = [a_{11}, a_{21}, a_{31}, a_{12}, a_{22}, a_{32}, a_{13}, a_{23}, a_{33}]

Suggested Fix

Standardized all matrices to use column-major indexing, as is standard in Julia/MATLAB convention:

# Column-major indexing: a[k] = a_{i,j} where k = i + (j-1)*n
A_row = zeros(n, n^2)
for i in 1:n
    for j in 1:n
        idx = i + (j-1)*n
        A_row[i, idx] = 1.0
    end
end

The sorting matrix T was also updated to use column-major and now has dimension (n-1) × n² instead of n × n²:

T = zeros(n-1, n^2)
for k in 1:n-1
    for j in 1:n
        idx_k = k + (j-1)*n
        idx_k1 = (k+1) + (j-1)*n
        T[k, idx_k] = 1.0
        T[k, idx_k1] = -1.0
    end
end

Issue 2: McCormick Envelope RHS Errors

Reference (Overleaf main.tex equations 19-22)

The McCormick envelopes for the bilinear term u_{ij} = a_{ij} · x_j where a_{ij} ∈ {0,1} and x_j ∈ [0, P_j]:

Eq. Constraint Description
(19) u_{ij} ≥ 0 Lower bound from a=0
(20) u_{ij} ≥ x̄_j · a_{ij} + x_j - x̄_j Lower bound from a=1
(21) u_{ij} ≤ x̄_j · a_{ij} Upper bound from x=x̄
(22) u_{ij} ≤ x_j Upper bound from a=0

Converting to F·y ≤ g·σ form

Envelope F row Correct RHS
(19) -I · u ≤ 0 zeros(n²)
(20) -I · u + A_P · a + x_j · x ≤ P P_out
(21) I · u - A_P · a ≤ 0 zeros(n²)
(22) I · u - x_j · x ≤ 0 zeros(n²)

Original Code

g = [ones(n), -ones(n), ones(n), -ones(n), zeros(n),
     zeros(n^2),   # envelope 1: correct
     -P_out,       # envelope 3: seems wrong (should be zeros?)
     zeros(n^2),   # envelope 2: seems wrong (should be P_out?)
     zeros(n^2)]   # envelope 4: correct

Fixed Code

g = [ones(n), -ones(n), ones(n), -ones(n), zeros(n-1),
     zeros(n^2),   # envelope 1: -u ≤ 0
     zeros(n^2),   # envelope 3: u - a·P ≤ 0 (Proposed Fix)
     P_out,        # envelope 2: -u + a·P + x ≤ P (Proposed Fix)
     zeros(n^2)]   # envelope 4: u - x ≤ 0

Issue 3: A_long Variable Reference

Problem

A_long = [A zeros(n,n^2) zeros(n,2*n)]  # A was undefined/row-major

Fix

A_long = [A_row zeros(n,n^2) zeros(n,2*n)]  # Uses column-major A_row

Files Modified

  1. src/implementation/palma_relaxation.jl

    • lin_palma() function (lines 14-145)
    • lin_palma_w_grad_input() function (lines 160-348)
  2. script/relaxed_palma_ratio.jl

    • Test script (lines 312-470)

Dimension Changes

Due to the T matrix correction, the constraint system dimensions changed:

Component Original Fixed
T matrix n × n² (n-1) × n²
g vector sorting block zeros(n) zeros(n-1)
Total F rows 4n + n + 4n² 4n + (n-1) + 4n²

Added Verification

Dimension assertions were added to catch future dimension mismatches:

@assert size(F, 2) == length(y) "F columns must match y length"
@assert size(F, 1) == length(g) "F rows must match g length"

Trying to align with overleaf formulation

The corrected implementation now matches the Overleaf formulation (assuming it is still accurate):

  1. Equations (5)-(6): Doubly stochastic constraints A·1 = 1, Aᵀ·1 = 1 are correctly implemented via A_row and A_col

  2. Equation (8): Sorting constraint S_i ≤ S_{i+1} is implemented via T · x_hat ≤ 0

  3. Equations (19)-(22): McCormick envelopes now have correct RHS values

  4. Equations (31)-(34): Charnes-Cooper transformation F·y ≤ g·σ with the currently written constraint structure


Testing

I created a proof-of-concept test script test/test_matrix_f.jl to verify:

  1. Matrix dimensions are consistent
  2. Permutation constraints are satisfied
  3. Sorting produces ascending order
  4. McCormick envelopes are correctly bounded
  5. Full constraint system F·y ≤ g is satisfied

Test Results

============================================================
Testing Matrix F Construction for Palma Ratio Linearization
============================================================

Test parameters: n = 3, P = [1.0, 2.0, 3.0]

--- Test 1: Dimension Consistency ---
Test Summary:         | Pass  Total
Dimension Consistency |    8      8

--- Test 2: Permutation Constraints ---
Identity permutation: Row sums = [1,1,1], Column sums = [1,1,1] ✓
Swap permutation: Row sums = [1,1,1], Column sums = [1,1,1] ✓
Test Summary:           | Pass  Total
Permutation Constraints |    2      2

--- Test 3: Sorting Constraint ---
Unsorted x: [3.0, 1.0, 2.0] → Sorted: [1.0, 2.0, 3.0]
T * x_hat = [-1.0, -1.0] (≤ 0) ✓
Test Summary:      | Pass  Total
Sorting Constraint |    2      2

--- Test 4: McCormick Envelope Bounds ---
For binary a ∈ {0,1}, bounds are tight at true value u = a·x ✓
Test Summary:       | Pass  Total
McCormick Envelopes |   12     12

--- Test 5: Full Constraint System ---
Max constraint violation: 0.0 ✓
All McCormick envelopes satisfied ✓
Test Summary:          | Pass  Total
Full Constraint System |   10     10

============================================================
All tests passed!
============================================================

Run the test with:

julia --project=. test/test_matrix_f.jl

Copy link

Copilot AI left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Pull request overview

This PR addresses three issues in the Palma ratio linearization implementation: standardizing vectorization indexing to column-major convention, correcting McCormick envelope RHS values, and fixing undefined variable references. The changes align the implementation with the mathematical formulation documented in Overleaf.

  • Standardized all matrix indexing to column-major convention (k = i + (j-1)*n) for consistency
  • Fixed McCormick envelope RHS values in the constraint vector g (swapped envelope 2 and 3 values)
  • Corrected sorting constraint matrix T to have (n-1) rows instead of n rows

Reviewed changes

Copilot reviewed 3 out of 3 changed files in this pull request and generated 7 comments.

File Description
test/test_matrix_f.jl New comprehensive test file verifying matrix dimensions, permutation constraints, sorting behavior, McCormick envelopes, and full constraint system
src/implementation/palma_relaxation.jl Updated both lin_palma() and lin_palma_w_grad_input() functions with column-major indexing, corrected T matrix dimensions, fixed McCormick RHS values, and added dimension assertions
script/relaxed_palma_ratio.jl Applied same corrections as implementation file to maintain consistency in the test script

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

@samtalki
Copy link
Collaborator Author

samtalki commented Jan 5, 2026

From the @copilot review, it looks like there was an inconsistency in -T vs T

  • lin_palma() used T, but lin_palma_w_grad_input() and relaxed_palma_ratio.jl used -T (incorrect).
  • In Overleaf T computes sorted[k] - sorted[k+1]. For ascending order (from eq 8: S_i ≤ S_{i+1}), we need T * x_hat ≤ 0, whereas -T would enforce descending order. I switched to ascending and used T*x_hat <= 0, let me know if this isn't what you wanted
  • Changed -T to T in two locations.

@copilot does this resolve your concern?

Copy link

Copilot AI commented Jan 5, 2026

@samtalki I've opened a new pull request, #5, to work on those changes. Once the pull request is ready, I'll request review from you.

samtalki and others added 6 commits January 14, 2026 14:18
Key changes:
- Add script/reformulation/load_shed_as_parameter.jl with new formulation
  where P_shed is a JuMP @expression (not @variable), eliminating
  equality constraints and simplifying the model structure
- Add script/reformulation/validate_reformulation.jl with tests
- Remove duplicate fairness functions from src/core/objective.jl
  (gini_index, jains_index, palma_ratio, alpha_fairness) that were
  causing precompilation errors due to method overwriting

The reformulation preserves full dynamic sorting via permutation matrix
optimization while using McCormick envelopes for bilinear terms and
Charnes-Cooper transformation for the ratio objective.

Co-Authored-By: Claude Opus 4.5 <noreply@anthropic.com>
DiffOpt fix:
- Increase regularization to 0.1 in MLD objective to keep pd variables
  interior, fixing sensitivity computation through KKT conditions
- Add regularization parameter to objective_fairly_weighted_max_load_served

Palma optimization fix:
- Add Taylor expansion constraint linking pshed_new to weight changes
- Fix objective from `Max, 0` to `Min, numerator_expr` (was maximizing zero!)
- Add Charnes-Cooper denominator normalization constraint
- Add trust region (±0.5) on weight changes for convergence stability
- Relax y_pshed lower bound to allow Taylor expansion flexibility
- Add MOI import and infeasibility handling
- Fix return statement to return optimized weights

Add diagnostic scripts:
- script/reformulation/debug/ with Jacobian analysis tools
- script/reformulation/test_bilevel_convergence.jl validates full pipeline

Verified: Palma ratio decreases from 5.97 to 5.16 (13.6% improvement)
over 5 bilevel iterations on ieee_13_aw test case.

Co-Authored-By: Claude Opus 4.5 <noreply@anthropic.com>
The McCormick relaxation (relax_binary=true) produces degenerate solutions
where the solver exploits loose envelopes to find objective=0. Binary
permutation is required for correct results.

Validated on IEEE 13-bus: NEW achieves 4.9% Palma improvement in 0.2s
vs OLD's 3.8% in 3.6s.

TODO: Investigate tighter McCormick cuts or alternative relaxations.

Co-Authored-By: Claude Opus 4.5 <noreply@anthropic.com>
Copy link
Owner

@awest32 awest32 left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The functionality seems to be here. Thank you so much for doing this!

@awest32 awest32 merged commit 368e716 into main Jan 16, 2026
2 of 3 checks passed
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants