use std::fmt::Write as FmtWrite;
use std::path::Path;
use surge_network::Network;
use surge_network::market::CostCurve;
use surge_network::network::BusType;
use thiserror::Error;
#[derive(Error, Debug)]
pub enum MatpowerWriteError {
#[error("I/O error: {0}")]
Io(#[from] std::io::Error),
#[error("format error: {0}")]
Fmt(#[from] std::fmt::Error),
}
pub fn write_file(network: &Network, path: &Path) -> Result<(), MatpowerWriteError> {
let content = to_string(network)?;
std::fs::write(path, content)?;
Ok(())
}
pub fn to_string(network: &Network) -> Result<String, MatpowerWriteError> {
let mut out = String::with_capacity(64 * 1024);
let fn_name = sanitize_name(&network.name);
writeln!(out, "function mpc = {fn_name}")?;
writeln!(
out,
"% Exported by Surge (https://github.com/amptimal/surge)"
)?;
writeln!(out, "mpc.version = '2';")?;
writeln!(out, "mpc.baseMVA = {};", network.base_mva)?;
writeln!(out)?;
let bus_map = network.bus_index_map();
let mut bus_demand_p = vec![0.0; network.buses.len()];
let mut bus_demand_q = vec![0.0; network.buses.len()];
for load in &network.loads {
if load.in_service {
if let Some(&idx) = bus_map.get(&load.bus) {
bus_demand_p[idx] += load.active_power_demand_mw;
bus_demand_q[idx] += load.reactive_power_demand_mvar;
}
}
}
writeln!(out, "mpc.bus = [")?;
writeln!(
out,
"%\tbus_i\ttype\tPd\tQd\tGs\tBs\tarea\tVm\tVa\tbaseKV\tzone\tVmax\tVmin"
)?;
for (bi, bus) in network.buses.iter().enumerate() {
let bus_type_code = match bus.bus_type {
BusType::PQ => 1,
BusType::PV => 2,
BusType::Slack => 3,
BusType::Isolated => 4,
};
let va_deg = bus.voltage_angle_rad.to_degrees();
writeln!(
out,
"\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{};",
bus.number,
bus_type_code,
bus_demand_p.get(bi).copied().unwrap_or(0.0),
bus_demand_q.get(bi).copied().unwrap_or(0.0),
bus.shunt_conductance_mw,
bus.shunt_susceptance_mvar,
bus.area,
bus.voltage_magnitude_pu,
va_deg,
bus.base_kv,
bus.zone,
bus.voltage_max_pu,
bus.voltage_min_pu
)?;
}
writeln!(out, "];")?;
writeln!(out)?;
writeln!(out, "mpc.gen = [")?;
writeln!(
out,
"%\tbus\tPg\tQg\tQmax\tQmin\tVg\tmBase\tstatus\tPmax\tPmin"
)?;
for g in &network.generators {
let status = if g.in_service { 1 } else { 0 };
let qmax = clamp_for_matpower(g.qmax, 9999.0);
let qmin = clamp_for_matpower(g.qmin, -9999.0);
let pmax = clamp_for_matpower(g.pmax, 9999.0);
let pmin = clamp_for_matpower(g.pmin, -9999.0);
let mbase = clamp_for_matpower(g.machine_base_mva, 100.0);
writeln!(
out,
"\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{};",
g.bus, g.p, g.q, qmax, qmin, g.voltage_setpoint_pu, mbase, status, pmax, pmin
)?;
}
for inj in &network.power_injections {
if !inj.in_service {
continue;
}
writeln!(
out,
"\t{}\t{}\t{}\t{}\t{}\t1.0\t100.0\t1\t{}\t{};",
inj.bus,
inj.active_power_injection_mw,
inj.reactive_power_injection_mvar,
inj.reactive_power_injection_mvar, inj.reactive_power_injection_mvar, inj.active_power_injection_mw, inj.active_power_injection_mw, )?;
}
writeln!(out, "];")?;
writeln!(out)?;
writeln!(out, "mpc.branch = [")?;
writeln!(
out,
"%\tfbus\ttbus\tr\tx\tb\trateA\trateB\trateC\ttap\tshift\tstatus\tangmin\tangmax"
)?;
for br in &network.branches {
let status = if br.in_service { 1 } else { 0 };
let tap_out = if (br.tap - 1.0).abs() < 1e-9 {
0.0
} else {
br.tap
};
let angmin = br
.angle_diff_min_rad
.unwrap_or(-2.0 * std::f64::consts::PI)
.to_degrees();
let angmax = br
.angle_diff_max_rad
.unwrap_or(2.0 * std::f64::consts::PI)
.to_degrees();
let ra = clamp_for_matpower(br.rating_a_mva, 0.0);
let rb = clamp_for_matpower(br.rating_b_mva, 0.0);
let rc = clamp_for_matpower(br.rating_c_mva, 0.0);
writeln!(
out,
"\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{};",
br.from_bus,
br.to_bus,
br.r,
br.x,
br.b,
ra,
rb,
rc,
tap_out,
br.phase_shift_rad.to_degrees(),
status,
angmin,
angmax
)?;
}
writeln!(out, "];")?;
writeln!(out)?;
let has_cost = network.generators.iter().any(|g| g.cost.is_some());
if has_cost {
writeln!(out, "mpc.gencost = [")?;
writeln!(out, "%\tmodel\tstartup\tshutdown\tn\tcost parameters")?;
for g in &network.generators {
match &g.cost {
Some(CostCurve::Polynomial {
startup,
shutdown,
coeffs,
}) => {
let n = coeffs.len();
write!(out, "\t2\t{startup}\t{shutdown}\t{n}")?;
for c in coeffs {
write!(out, "\t{c}")?;
}
writeln!(out, ";")?;
}
Some(CostCurve::PiecewiseLinear {
startup,
shutdown,
points,
}) => {
let n = points.len();
write!(out, "\t1\t{startup}\t{shutdown}\t{n}")?;
for (p, c) in points {
write!(out, "\t{p}\t{c}")?;
}
writeln!(out, ";")?;
}
None => {
writeln!(out, "\t2\t0\t0\t2\t0\t0;")?;
}
}
}
for inj in &network.power_injections {
if !inj.in_service {
continue;
}
writeln!(out, "\t2\t0\t0\t1\t0;")?;
}
writeln!(out, "];")?;
}
if network.hvdc.has_explicit_dc_topology() {
writeln!(out)?;
writeln!(out, "mpc.busdc = [")?;
writeln!(
out,
"%\tbusdc_i\tgrid\tPdc\tVdc\tbasekVdc\tVdcmax\tVdcmin\tCdc"
)?;
for dc_grid in &network.hvdc.dc_grids {
for b in &dc_grid.buses {
writeln!(
out,
"\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{};",
b.bus_id,
dc_grid.id,
b.p_dc_mw,
b.v_dc_pu,
b.base_kv_dc,
b.v_dc_max,
b.v_dc_min,
b.cost
)?;
}
}
writeln!(out, "];")?;
writeln!(out)?;
writeln!(out, "mpc.convdc = [")?;
writeln!(
out,
"%\tdc_bus\tac_bus\ttype_dc\ttype_ac\tP_g\tQ_g\tislcc\tVtar\trtf\txtf\ttransformer\ttm\tbf\tfilter\trc\txc\treactor\tbasekVac\tVmmax\tVmmin\tImax\tstatus\tLossA\tLossB\tLossCrec\tLossCinv\tdroop\tPdcset\tVdcset\tPacmax\tPacmin\tQacmax\tQacmin"
)?;
for dc_grid in &network.hvdc.dc_grids {
for converter in &dc_grid.converters {
let Some(c) = converter.as_vsc() else {
continue;
};
let is_lcc = if c.is_lcc { 1 } else { 0 };
let transformer = if c.transformer { 1 } else { 0 };
let filter = if c.filter { 1 } else { 0 };
let reactor = if c.reactor { 1 } else { 0 };
let status = if c.status { 1 } else { 0 };
let pac_max = clamp_for_matpower(c.active_power_ac_max_mw, 9999.0);
let pac_min = clamp_for_matpower(c.active_power_ac_min_mw, -9999.0);
let qac_max = clamp_for_matpower(c.reactive_power_ac_max_mvar, 9999.0);
let qac_min = clamp_for_matpower(c.reactive_power_ac_min_mvar, -9999.0);
writeln!(
out,
"\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{};",
c.dc_bus,
c.ac_bus,
c.control_type_dc,
c.control_type_ac,
c.active_power_mw,
c.reactive_power_mvar,
is_lcc,
c.voltage_setpoint_pu,
c.transformer_r_pu,
c.transformer_x_pu,
transformer,
c.tap_ratio,
c.filter_susceptance_pu,
filter,
c.reactor_r_pu,
c.reactor_x_pu,
reactor,
c.base_kv_ac,
c.voltage_max_pu,
c.voltage_min_pu,
c.current_max_pu,
status,
c.loss_constant_mw,
c.loss_linear,
c.loss_quadratic_rectifier,
c.loss_quadratic_inverter,
c.droop,
c.power_dc_setpoint_mw,
c.voltage_dc_setpoint_pu,
pac_max,
pac_min,
qac_max,
qac_min
)?;
}
}
writeln!(out, "];")?;
writeln!(out)?;
writeln!(out, "mpc.branchdc = [")?;
writeln!(
out,
"%\tfbusdc\ttbusdc\tr\tl\tc\trateA\trateB\trateC\tstatus"
)?;
for dc_grid in &network.hvdc.dc_grids {
for br in &dc_grid.branches {
let status = if br.status { 1 } else { 0 };
writeln!(
out,
"\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{};",
br.from_bus,
br.to_bus,
br.r_ohm,
br.l_mh,
br.c_uf,
br.rating_a_mva,
br.rating_b_mva,
br.rating_c_mva,
status
)?;
}
}
writeln!(out, "];")?;
}
Ok(out)
}
fn clamp_for_matpower(v: f64, fallback: f64) -> f64 {
if !v.is_finite() || v >= f64::MAX / 2.0 || v <= f64::MIN / 2.0 {
fallback
} else {
v
}
}
fn sanitize_name(name: &str) -> String {
let s: String = name
.chars()
.map(|c| {
if c.is_alphanumeric() || c == '_' {
c
} else {
'_'
}
})
.collect();
if s.starts_with(|c: char| c.is_ascii_digit()) {
format!("case_{s}")
} else if s.is_empty() {
"case_unknown".to_string()
} else {
s
}
}
#[cfg(test)]
mod tests {
use super::*;
use surge_network::Network;
use surge_network::market::CostCurve;
use surge_network::network::{
Branch, Bus, BusType, DcBranch, DcBus, DcConverterStation, Generator, Load,
};
fn simple_network() -> Network {
let mut net = Network::new("test_case");
net.base_mva = 100.0;
let mut slack = Bus::new(1, BusType::Slack, 345.0);
slack.voltage_magnitude_pu = 1.04;
slack.voltage_angle_rad = 0.0;
net.buses.push(slack);
let pq = Bus::new(2, BusType::PQ, 345.0);
net.buses.push(pq);
net.loads.push(Load::new(2, 100.0, 35.0));
let mut g = Generator::new(1, 80.0, 1.04);
g.qmax = 300.0;
g.qmin = -300.0;
g.pmax = 250.0;
g.pmin = 10.0;
g.cost = Some(CostCurve::Polynomial {
startup: 1500.0,
shutdown: 0.0,
coeffs: vec![0.11, 5.0, 150.0],
});
net.generators.push(g);
net.branches
.push(Branch::new_line(1, 2, 0.01938, 0.05917, 0.0528));
net
}
#[test]
fn test_write_produces_matpower_sections() {
let net = simple_network();
let s = to_string(&net).unwrap();
assert!(s.contains("function mpc = test_case"));
assert!(s.contains("mpc.baseMVA = 100"));
assert!(s.contains("mpc.bus = ["));
assert!(s.contains("mpc.gen = ["));
assert!(s.contains("mpc.branch = ["));
assert!(s.contains("mpc.gencost = ["));
}
#[test]
fn test_roundtrip_matpower() {
use crate::matpower::parse_str;
let net = simple_network();
let s = to_string(&net).unwrap();
let net2 = parse_str(&s).expect("round-trip parse failed");
assert_eq!(net2.n_buses(), net.n_buses());
assert_eq!(net2.n_branches(), net.n_branches());
assert!((net2.base_mva - net.base_mva).abs() < 1e-6);
assert!(
(net2.buses[0].voltage_magnitude_pu - net.buses[0].voltage_magnitude_pu).abs() < 1e-6
);
let pd1 = net.bus_load_p_mw();
let pd2 = net2.bus_load_p_mw();
assert!((pd2[1] - pd1[1]).abs() < 1e-6);
}
#[test]
fn test_tap_convention() {
let net = simple_network();
let s = to_string(&net).unwrap();
assert!(s.contains("mpc.branch"));
}
#[test]
fn test_va_in_degrees() {
let mut net = Network::new("va_test");
net.base_mva = 100.0;
let mut bus = Bus::new(1, BusType::Slack, 138.0);
bus.voltage_angle_rad = std::f64::consts::PI / 6.0; net.buses.push(bus);
net.generators.push(Generator::new(1, 0.0, 1.0));
let s = to_string(&net).unwrap();
assert!(
!s.contains("0.5235"),
"va written in radians, expected degrees"
);
let va_deg_str = format!("{}", (std::f64::consts::PI / 6.0_f64).to_degrees());
assert!(
s.contains(&va_deg_str[..4]), "expected va in degrees (~30) in output, got: {s}"
);
}
#[test]
fn test_piecewise_linear_gencost() {
let mut net = Network::new("pwl_case");
net.base_mva = 100.0;
net.buses.push(Bus::new(1, BusType::Slack, 138.0));
let mut g = Generator::new(1, 0.0, 1.0);
g.cost = Some(CostCurve::PiecewiseLinear {
startup: 0.0,
shutdown: 0.0,
points: vec![(0.0, 0.0), (100.0, 1000.0), (200.0, 3000.0)],
});
net.generators.push(g);
let s = to_string(&net).unwrap();
assert!(s.contains("\t1\t"));
}
#[test]
fn test_file_roundtrip() {
let net = simple_network();
let tmp = std::env::temp_dir().join("surge_mw_writer_test.m");
write_file(&net, &tmp).unwrap();
let content = std::fs::read_to_string(&tmp).unwrap();
assert!(content.contains("mpc.bus = ["));
let _ = std::fs::remove_file(&tmp);
}
#[test]
fn test_power_injection_rows_get_matching_gencost_entries() {
let mut net = simple_network();
net.power_injections
.push(surge_network::network::PowerInjection::new(2, 15.0, 3.0));
let s = to_string(&net).unwrap();
let gen_section = s
.split("mpc.gen = [")
.nth(1)
.and_then(|part| part.split("];").next())
.unwrap();
let gencost_section = s
.split("mpc.gencost = [")
.nth(1)
.and_then(|part| part.split("];").next())
.unwrap();
let gen_count = gen_section
.lines()
.filter(|line| line.starts_with('\t') && line.trim_end().ends_with(';'))
.count();
let gencost_count = gencost_section
.lines()
.filter(|line| line.starts_with('\t') && line.trim_end().ends_with(';'))
.count();
assert_eq!(gen_count, 2);
assert_eq!(gencost_count, 2);
}
fn dc_network() -> Network {
let mut net = simple_network();
let grid = net.hvdc.ensure_dc_grid(1, None);
grid.buses.push(DcBus {
bus_id: 1,
p_dc_mw: 0.0,
v_dc_pu: 1.0,
base_kv_dc: 400.0,
v_dc_max: 1.1,
v_dc_min: 0.9,
cost: 0.0,
g_shunt_siemens: 0.0,
r_ground_ohm: 0.0,
});
grid.buses.push(DcBus {
bus_id: 2,
p_dc_mw: 0.0,
v_dc_pu: 1.0,
base_kv_dc: 400.0,
v_dc_max: 1.1,
v_dc_min: 0.9,
cost: 0.0,
g_shunt_siemens: 0.0,
r_ground_ohm: 0.0,
});
grid.converters.push(
DcConverterStation {
id: String::new(),
dc_bus: 1,
ac_bus: 1,
control_type_dc: 2,
control_type_ac: 1,
active_power_mw: 200.0,
reactive_power_mvar: 30.0,
is_lcc: false,
voltage_setpoint_pu: 1.0,
transformer_r_pu: 0.0015,
transformer_x_pu: 0.1121,
transformer: true,
tap_ratio: 1.0,
filter_susceptance_pu: 0.0887,
filter: true,
reactor_r_pu: 0.0001,
reactor_x_pu: 0.16428,
reactor: true,
base_kv_ac: 345.0,
voltage_max_pu: 1.1,
voltage_min_pu: 0.9,
current_max_pu: 1.1,
status: true,
loss_constant_mw: 1.103,
loss_linear: 0.887,
loss_quadratic_rectifier: 2.885,
loss_quadratic_inverter: 2.885,
droop: 0.0,
power_dc_setpoint_mw: 200.0,
voltage_dc_setpoint_pu: 1.0,
active_power_ac_max_mw: 500.0,
active_power_ac_min_mw: -500.0,
reactive_power_ac_max_mvar: 300.0,
reactive_power_ac_min_mvar: -300.0,
}
.into(),
);
grid.converters.push(
DcConverterStation {
id: String::new(),
dc_bus: 2,
ac_bus: 2,
control_type_dc: 1,
control_type_ac: 1,
active_power_mw: -195.0,
reactive_power_mvar: 20.0,
is_lcc: false,
voltage_setpoint_pu: 1.0,
transformer_r_pu: 0.0015,
transformer_x_pu: 0.1121,
transformer: true,
tap_ratio: 1.0,
filter_susceptance_pu: 0.0887,
filter: true,
reactor_r_pu: 0.0001,
reactor_x_pu: 0.16428,
reactor: true,
base_kv_ac: 345.0,
voltage_max_pu: 1.1,
voltage_min_pu: 0.9,
current_max_pu: 1.1,
status: true,
loss_constant_mw: 1.103,
loss_linear: 0.887,
loss_quadratic_rectifier: 2.885,
loss_quadratic_inverter: 2.885,
droop: 0.0,
power_dc_setpoint_mw: -195.0,
voltage_dc_setpoint_pu: 1.0,
active_power_ac_max_mw: 500.0,
active_power_ac_min_mw: -500.0,
reactive_power_ac_max_mvar: 300.0,
reactive_power_ac_min_mvar: -300.0,
}
.into(),
);
grid.branches.push(DcBranch {
id: "dc_branch_1".into(),
from_bus: 1,
to_bus: 2,
r_ohm: 5.0,
l_mh: 25.0,
c_uf: 12.0,
rating_a_mva: 500.0,
rating_b_mva: 0.0,
rating_c_mva: 0.0,
status: true,
});
net
}
#[test]
fn test_dc_sections_written() {
let net = dc_network();
let s = to_string(&net).unwrap();
assert!(s.contains("mpc.busdc = ["));
assert!(s.contains("mpc.convdc = ["));
assert!(s.contains("mpc.branchdc = ["));
}
#[test]
fn test_dc_roundtrip() {
use crate::matpower::parse_str;
let net = dc_network();
let s = to_string(&net).unwrap();
let net2 = parse_str(&s).expect("DC round-trip parse failed");
let dc_buses: Vec<_> = net2.hvdc.dc_buses().collect();
let dc_converters: Vec<_> = net2
.hvdc
.dc_converters()
.filter_map(|c| c.as_vsc())
.collect();
let dc_branches: Vec<_> = net2.hvdc.dc_branches().collect();
assert_eq!(dc_buses.len(), 2);
assert_eq!(dc_buses[0].bus_id, 1);
assert!((dc_buses[0].base_kv_dc - 400.0).abs() < 1e-6);
assert_eq!(dc_converters.len(), 2);
assert_eq!(dc_converters[0].control_type_dc, 2);
assert!((dc_converters[0].active_power_mw - 200.0).abs() < 1e-6);
assert!((dc_converters[0].loss_constant_mw - 1.103).abs() < 1e-6);
assert!((dc_converters[0].loss_linear - 0.887).abs() < 1e-6);
assert!((dc_converters[0].loss_quadratic_rectifier - 2.885).abs() < 1e-6);
assert!((dc_converters[0].active_power_ac_max_mw - 500.0).abs() < 1e-6);
assert!((dc_converters[0].active_power_ac_min_mw - (-500.0)).abs() < 1e-6);
assert!((dc_converters[0].voltage_dc_setpoint_pu - 1.0).abs() < 1e-6);
assert!(dc_converters[0].transformer);
assert!(dc_converters[0].filter);
assert!(dc_converters[0].reactor);
assert_eq!(dc_branches.len(), 1);
assert!((dc_branches[0].r_ohm - 5.0).abs() < 1e-6);
assert!((dc_branches[0].l_mh - 25.0).abs() < 1e-6);
assert!((dc_branches[0].c_uf - 12.0).abs() < 1e-6);
assert!(dc_branches[0].status);
}
#[test]
fn test_no_dc_sections_when_empty() {
let net = simple_network();
let s = to_string(&net).unwrap();
assert!(!s.contains("mpc.busdc"));
assert!(!s.contains("mpc.convdc"));
assert!(!s.contains("mpc.branchdc"));
}
}