use std::collections::HashMap;
use cobre_core::entities::hydro::HydroGenerationModel;
use cobre_core::{
Bus, CascadeTopology, ConstraintSense, EntityId, GenericConstraint, Hydro, Line, LoadModel,
NonControllableSource, ResolvedBounds, ResolvedExchangeFactors,
ResolvedGenericConstraintBounds, ResolvedLoadFactors, ResolvedNcsBounds, ResolvedNcsFactors,
ResolvedPenalties, Stage, System, Thermal,
};
use crate::generic_constraints::resolve_variable_ref;
use cobre_solver::StageTemplate;
use cobre_stochastic::normal::precompute::PrecomputedNormal;
use cobre_stochastic::par::precompute::PrecomputedPar;
use crate::error::SddpError;
use crate::hydro_models::{
EvaporationModel, EvaporationModelSet, ProductionModelSet, ResolvedProductionModel,
};
use crate::indexer::StageIndexer;
use crate::inflow_method::InflowNonNegativityMethod;
#[derive(Debug, Clone)]
pub struct PatchBuffer {
pub indices: Vec<usize>,
pub lower: Vec<f64>,
pub upper: Vec<f64>,
hydro_count: usize,
max_par_order: usize,
load_bus_count: usize,
max_blocks: usize,
active_load_patches: usize,
active_z_inflow_patches: usize,
}
impl PatchBuffer {
#[must_use]
pub fn new(
hydro_count: usize,
max_par_order: usize,
n_load_buses: usize,
max_blocks: usize,
) -> Self {
let capacity = hydro_count * (2 + max_par_order) + n_load_buses * max_blocks + hydro_count;
Self {
indices: vec![0; capacity],
lower: vec![0.0; capacity],
upper: vec![0.0; capacity],
hydro_count,
max_par_order,
load_bus_count: n_load_buses,
max_blocks,
active_load_patches: 0,
active_z_inflow_patches: 0,
}
}
pub fn fill_forward_patches(
&mut self,
indexer: &StageIndexer,
state: &[f64],
noise: &[f64],
base_row: usize,
row_scale: &[f64],
) {
debug_assert_eq!(
state.len(),
indexer.n_state,
"state slice length {got} != n_state {expected}",
got = state.len(),
expected = indexer.n_state,
);
debug_assert!(
noise.len() == indexer.hydro_count || noise.is_empty(),
"noise slice length {got} must equal hydro_count {expected} or be empty",
got = noise.len(),
expected = indexer.hydro_count,
);
let n = self.hydro_count;
let l = self.max_par_order;
for (h, &sv) in state[..n].iter().enumerate() {
self.indices[h] = h;
let scaled = if row_scale.is_empty() {
sv
} else {
sv * row_scale[h]
};
self.lower[h] = scaled;
self.upper[h] = scaled;
}
for lag in 0..l {
for h in 0..n {
let slot = n + lag * n + h;
self.indices[slot] = slot;
let sv = state[slot];
let scaled = if row_scale.is_empty() {
sv
} else {
sv * row_scale[slot]
};
self.lower[slot] = scaled;
self.upper[slot] = scaled;
}
}
let cat3_start = n * (1 + l);
for (h, &nv) in noise.iter().enumerate() {
let slot = cat3_start + h;
self.indices[slot] = ar_dynamics_row_offset(base_row, h);
self.lower[slot] = nv;
self.upper[slot] = nv;
}
}
pub fn fill_state_patches(&mut self, indexer: &StageIndexer, state: &[f64], row_scale: &[f64]) {
debug_assert_eq!(
state.len(),
indexer.n_state,
"state slice length {got} != n_state {expected}",
got = state.len(),
expected = indexer.n_state,
);
let n = self.hydro_count;
let l = self.max_par_order;
for (h, &sv) in state[..n].iter().enumerate() {
self.indices[h] = h;
let scaled = if row_scale.is_empty() {
sv
} else {
sv * row_scale[h]
};
self.lower[h] = scaled;
self.upper[h] = scaled;
}
for lag in 0..l {
for h in 0..n {
let slot = n + lag * n + h;
self.indices[slot] = slot;
let sv = state[slot];
let scaled = if row_scale.is_empty() {
sv
} else {
sv * row_scale[slot]
};
self.lower[slot] = scaled;
self.upper[slot] = scaled;
}
}
}
pub fn fill_load_patches(
&mut self,
load_row_start: usize,
n_blocks: usize,
load_rhs: &[f64],
bus_positions: &[usize],
row_scale: &[f64],
) {
debug_assert_eq!(
load_rhs.len(),
self.load_bus_count * n_blocks,
"load_rhs length {got} != load_bus_count*n_blocks {expected}",
got = load_rhs.len(),
expected = self.load_bus_count * n_blocks,
);
debug_assert_eq!(
bus_positions.len(),
self.load_bus_count,
"bus_positions length {got} != load_bus_count {expected}",
got = bus_positions.len(),
expected = self.load_bus_count,
);
debug_assert!(
n_blocks <= self.max_blocks,
"n_blocks {n_blocks} exceeds max_blocks {mb}",
mb = self.max_blocks,
);
let cat4_start = self.hydro_count * (2 + self.max_par_order);
let mut slot = cat4_start;
for (i, &bus_pos) in bus_positions.iter().enumerate() {
for blk in 0..n_blocks {
let row = load_row_start + bus_pos * n_blocks + blk;
let rhs = load_rhs[i * n_blocks + blk];
let scaled = if row_scale.is_empty() {
rhs
} else {
rhs * row_scale[row]
};
self.indices[slot] = row;
self.lower[slot] = scaled;
self.upper[slot] = scaled;
slot += 1;
}
}
self.active_load_patches = self.load_bus_count * n_blocks;
}
pub fn fill_z_inflow_patches(
&mut self,
z_inflow_row_start: usize,
z_inflow_rhs: &[f64],
row_scale: &[f64],
) {
let n = self.hydro_count;
if n == 0 || z_inflow_rhs.is_empty() {
self.active_z_inflow_patches = 0;
return;
}
let cat5_start = self.hydro_count * (2 + self.max_par_order) + self.active_load_patches;
for (h, &rhs) in z_inflow_rhs.iter().enumerate().take(n) {
let slot = cat5_start + h;
let row = z_inflow_row_start + h;
let scaled = if row_scale.is_empty() {
rhs
} else {
rhs * row_scale[row]
};
self.indices[slot] = row;
self.lower[slot] = scaled;
self.upper[slot] = scaled;
}
self.active_z_inflow_patches = n;
}
#[must_use]
#[inline]
pub fn forward_patch_count(&self) -> usize {
self.hydro_count * (2 + self.max_par_order)
+ self.active_load_patches
+ self.active_z_inflow_patches
}
#[must_use]
#[inline]
pub fn state_patch_count(&self) -> usize {
self.hydro_count * (1 + self.max_par_order)
}
}
#[must_use]
#[inline]
pub fn ar_dynamics_row_offset(base_row: usize, hydro_index: usize) -> usize {
base_row + hydro_index
}
#[derive(Debug, Clone)]
pub struct StageTemplates {
pub templates: Vec<StageTemplate>,
pub base_rows: Vec<usize>,
pub noise_scale: Vec<f64>,
pub zeta_per_stage: Vec<f64>,
pub block_hours_per_stage: Vec<Vec<f64>>,
pub n_hydros: usize,
pub load_balance_row_starts: Vec<usize>,
pub n_load_buses: usize,
pub load_bus_indices: Vec<usize>,
pub generic_constraint_row_entries: Vec<Vec<GenericConstraintRowEntry>>,
pub ncs_col_starts: Vec<usize>,
pub n_ncs_per_stage: Vec<usize>,
pub active_ncs_indices: Vec<Vec<usize>>,
pub diversion_upstream: HashMap<EntityId, Vec<usize>>,
pub hydro_productivities_per_stage: Vec<Vec<f64>>,
}
const M3S_TO_HM3: f64 = 3_600.0 / 1_000_000.0;
pub(crate) const COST_SCALE_FACTOR: f64 = 1_000.0;
pub(crate) const Q_EV_SAFETY_MARGIN: f64 = 2.0;
pub(crate) const OVER_EVAPORATION_COST_MULTIPLIER: f64 = 100.0;
#[derive(Debug, Clone)]
pub struct GenericConstraintRowEntry {
pub constraint_idx: usize,
pub entity_id: i32,
pub block_idx: usize,
pub bound: f64,
pub sense: ConstraintSense,
pub slack_enabled: bool,
pub slack_penalty: f64,
pub slack_plus_col: Option<usize>,
pub slack_minus_col: Option<usize>,
}
struct TemplateBuildCtx<'a> {
hydros: &'a [Hydro],
thermals: &'a [Thermal],
lines: &'a [Line],
buses: &'a [Bus],
load_models: &'a [LoadModel],
cascade: &'a CascadeTopology,
bounds: &'a ResolvedBounds,
penalties: &'a ResolvedPenalties,
hydro_pos: HashMap<EntityId, usize>,
thermal_pos: HashMap<EntityId, usize>,
line_pos: HashMap<EntityId, usize>,
bus_pos: HashMap<EntityId, usize>,
inflow_method: &'a InflowNonNegativityMethod,
par_lp: &'a PrecomputedPar,
production_models: &'a ProductionModelSet,
evaporation_models: &'a EvaporationModelSet,
generic_constraints: &'a [GenericConstraint],
resolved_generic_bounds: &'a ResolvedGenericConstraintBounds,
resolved_load_factors: &'a ResolvedLoadFactors,
resolved_exchange_factors: &'a ResolvedExchangeFactors,
non_controllable_sources: &'a [NonControllableSource],
resolved_ncs_bounds: &'a ResolvedNcsBounds,
resolved_ncs_factors: &'a ResolvedNcsFactors,
diversion_upstream: HashMap<EntityId, Vec<usize>>,
n_hydros: usize,
n_thermals: usize,
n_lines: usize,
n_buses: usize,
max_par_order: usize,
has_penalty: bool,
}
struct StageLayout {
n_blks: usize,
n_h: usize,
lag_order: usize,
col_turbine_start: usize,
col_spillage_start: usize,
col_diversion_start: usize,
col_thermal_start: usize,
col_line_fwd_start: usize,
col_line_rev_start: usize,
col_deficit_start: usize,
max_deficit_segments: usize,
col_excess_start: usize,
col_inflow_slack_start: usize,
col_generation_start: usize,
row_water_balance_start: usize,
row_load_balance_start: usize,
row_fpha_start: usize,
row_evap_start: usize,
col_evap_start: usize,
col_withdrawal_slack_start: usize,
col_ncs_start: usize,
n_ncs: usize,
active_ncs_indices: Vec<usize>,
num_cols: usize,
row_generic_start: usize,
num_rows: usize,
n_generic_rows: usize,
row_z_inflow_start: usize,
col_z_inflow_start: usize,
n_state: usize,
n_dual_relevant: usize,
zeta: f64,
fpha_hydro_indices: Vec<usize>,
fpha_planes_per_hydro: Vec<usize>,
evap_hydro_indices: Vec<usize>,
pub(crate) generic_constraint_rows: Vec<GenericConstraintRowEntry>,
}
impl StageLayout {
#[allow(clippy::too_many_lines)]
fn new(ctx: &TemplateBuildCtx<'_>, stage: &Stage, stage_idx: usize) -> Self {
let n_blks = stage.blocks.len();
let idx = StageIndexer::new(ctx.n_hydros, ctx.max_par_order);
let decision_start = idx.theta + 1;
let col_turbine_start = decision_start;
let col_spillage_start = col_turbine_start + ctx.n_hydros * n_blks;
let col_diversion_start = col_spillage_start + ctx.n_hydros * n_blks;
let col_thermal_start = col_diversion_start + ctx.n_hydros * n_blks;
let col_line_fwd_start = col_thermal_start + ctx.n_thermals * n_blks;
let col_line_rev_start = col_line_fwd_start + ctx.n_lines * n_blks;
let max_deficit_segments = ctx
.buses
.iter()
.map(|b| b.deficit_segments.len())
.max()
.unwrap_or(0);
let col_deficit_start = col_line_rev_start + ctx.n_lines * n_blks;
let col_excess_start = col_deficit_start + ctx.n_buses * max_deficit_segments * n_blks;
let col_excess_end = col_excess_start + ctx.n_buses * n_blks;
let col_inflow_slack_start = col_excess_end;
let n_slack_cols = if ctx.has_penalty { ctx.n_hydros } else { 0 };
let mut fpha_hydro_indices: Vec<usize> = Vec::new();
let mut fpha_planes_per_hydro: Vec<usize> = Vec::new();
for h_idx in 0..ctx.n_hydros {
if let ResolvedProductionModel::Fpha { planes, .. } =
ctx.production_models.model(h_idx, stage_idx)
{
fpha_hydro_indices.push(h_idx);
fpha_planes_per_hydro.push(planes.len());
}
}
let n_fpha_hydros = fpha_hydro_indices.len();
let col_generation_start = col_inflow_slack_start + n_slack_cols;
let n_generation_cols = n_fpha_hydros * n_blks;
let col_generation_end = col_generation_start + n_generation_cols;
let mut evap_hydro_indices: Vec<usize> = Vec::new();
for h_idx in 0..ctx.n_hydros {
if matches!(
ctx.evaporation_models.model(h_idx),
EvaporationModel::Linearized { .. }
) {
evap_hydro_indices.push(h_idx);
}
}
let n_evap_hydros = evap_hydro_indices.len();
let col_evap_start = col_generation_end;
let n_evap_cols = n_evap_hydros * 3;
let col_withdrawal_slack_start = col_evap_start + n_evap_cols;
let withdrawal_slack_end = col_withdrawal_slack_start + ctx.n_hydros;
let n_state = idx.n_state;
let n_dual_relevant = n_state;
let row_water_balance_start = n_dual_relevant + ctx.n_hydros;
let row_load_balance_start = row_water_balance_start + ctx.n_hydros;
let row_load_balance_end = row_load_balance_start + ctx.n_buses * n_blks;
let row_fpha_start = row_load_balance_end;
let n_fpha_rows: usize = fpha_planes_per_hydro.iter().map(|&p| p * n_blks).sum();
let row_evap_start = row_fpha_start + n_fpha_rows;
let evap_rows_end = row_evap_start + n_evap_hydros;
let mut active_ncs_indices: Vec<usize> = Vec::new();
for (ncs_idx, ncs) in ctx.non_controllable_sources.iter().enumerate() {
let entered = ncs.entry_stage_id.is_none_or(|entry| entry <= stage.id);
let not_exited = ncs.exit_stage_id.is_none_or(|exit| stage.id < exit);
if entered && not_exited {
active_ncs_indices.push(ncs_idx);
}
}
let n_active_ncs = active_ncs_indices.len();
let col_ncs_start = withdrawal_slack_end;
let n_ncs_cols = n_active_ncs * n_blks;
let col_ncs_end = col_ncs_start + n_ncs_cols;
let row_generic_start = evap_rows_end;
let col_generic_slack_start = col_ncs_end;
let mut n_generic_rows: usize = 0;
let mut n_generic_slack_cols: usize = 0;
let mut generic_constraint_rows: Vec<GenericConstraintRowEntry> = Vec::new();
for (constraint_idx, constraint) in ctx.generic_constraints.iter().enumerate() {
if !ctx
.resolved_generic_bounds
.is_active(constraint_idx, stage.id)
{
continue;
}
let bound_entries = ctx
.resolved_generic_bounds
.bounds_for_stage(constraint_idx, stage.id);
for &(block_id, bound) in bound_entries {
match block_id {
None => {
for block_idx in 0..n_blks {
let row_offset = row_generic_start + n_generic_rows;
let (slack_plus_col, slack_minus_col) = if constraint.slack.enabled {
let plus_col = col_generic_slack_start + n_generic_slack_cols;
n_generic_slack_cols += 1;
let minus_col = if constraint.sense == ConstraintSense::Equal {
let mc = col_generic_slack_start + n_generic_slack_cols;
n_generic_slack_cols += 1;
Some(mc)
} else {
None
};
(Some(plus_col), minus_col)
} else {
(None, None)
};
let _ = row_offset; n_generic_rows += 1;
generic_constraint_rows.push(GenericConstraintRowEntry {
constraint_idx,
entity_id: constraint.id.0,
block_idx,
bound,
sense: constraint.sense,
slack_enabled: constraint.slack.enabled,
slack_penalty: constraint.slack.penalty.unwrap_or(0.0),
slack_plus_col,
slack_minus_col,
});
}
}
Some(blk_id) => {
#[allow(clippy::cast_sign_loss)]
let block_idx = blk_id as usize;
let (slack_plus_col, slack_minus_col) = if constraint.slack.enabled {
let plus_col = col_generic_slack_start + n_generic_slack_cols;
n_generic_slack_cols += 1;
let minus_col = if constraint.sense == ConstraintSense::Equal {
let mc = col_generic_slack_start + n_generic_slack_cols;
n_generic_slack_cols += 1;
Some(mc)
} else {
None
};
(Some(plus_col), minus_col)
} else {
(None, None)
};
n_generic_rows += 1;
generic_constraint_rows.push(GenericConstraintRowEntry {
constraint_idx,
entity_id: constraint.id.0,
block_idx,
bound,
sense: constraint.sense,
slack_enabled: constraint.slack.enabled,
slack_penalty: constraint.slack.penalty.unwrap_or(0.0),
slack_plus_col,
slack_minus_col,
});
}
}
}
}
let col_z_inflow_start = idx.z_inflow.start; let row_z_inflow_start = idx.z_inflow_row_start;
let num_cols = col_generic_slack_start + n_generic_slack_cols;
let num_rows = row_generic_start + n_generic_rows;
let total_stage_hours: f64 = stage.blocks.iter().map(|b| b.duration_hours).sum();
let zeta = total_stage_hours * M3S_TO_HM3;
Self {
n_blks,
n_h: ctx.n_hydros,
lag_order: ctx.max_par_order,
col_turbine_start,
col_spillage_start,
col_diversion_start,
col_thermal_start,
col_line_fwd_start,
col_line_rev_start,
col_deficit_start,
max_deficit_segments,
col_excess_start,
col_inflow_slack_start,
col_generation_start,
col_evap_start,
col_withdrawal_slack_start,
col_ncs_start,
n_ncs: n_active_ncs,
active_ncs_indices,
num_cols,
row_water_balance_start,
row_load_balance_start,
row_fpha_start,
row_evap_start,
row_generic_start,
num_rows,
n_generic_rows,
row_z_inflow_start,
col_z_inflow_start,
n_state,
n_dual_relevant,
zeta,
fpha_hydro_indices,
fpha_planes_per_hydro,
evap_hydro_indices,
generic_constraint_rows,
}
}
}
#[allow(clippy::too_many_lines)]
fn fill_stage_columns(
ctx: &TemplateBuildCtx<'_>,
stage: &Stage,
stage_idx: usize,
layout: &StageLayout,
) -> (Vec<f64>, Vec<f64>, Vec<f64>) {
let idx = StageIndexer::new(ctx.n_hydros, ctx.max_par_order);
let mut col_lower = vec![0.0_f64; layout.num_cols];
let mut col_upper = vec![f64::INFINITY; layout.num_cols];
let mut objective = vec![0.0_f64; layout.num_cols];
for (h_idx, _hydro) in ctx.hydros.iter().enumerate() {
let hb = ctx.bounds.hydro_bounds(h_idx, stage_idx);
col_lower[h_idx] = hb.min_storage_hm3;
col_upper[h_idx] = hb.max_storage_hm3;
col_lower[idx.storage_in.start + h_idx] = f64::NEG_INFINITY;
col_upper[idx.storage_in.start + h_idx] = f64::INFINITY;
}
for lag_col in idx.inflow_lags.clone() {
col_lower[lag_col] = f64::NEG_INFINITY;
col_upper[lag_col] = f64::INFINITY;
}
col_lower[idx.theta] = 0.0;
col_upper[idx.theta] = f64::INFINITY;
objective[idx.theta] = 1.0;
for (h_idx, _hydro) in ctx.hydros.iter().enumerate() {
let hb = ctx.bounds.hydro_bounds(h_idx, stage_idx);
let hp = ctx.penalties.hydro_penalties(h_idx, stage_idx);
let model = ctx.production_models.model(h_idx, stage_idx);
let is_fpha = matches!(model, ResolvedProductionModel::Fpha { .. });
let turb_upper = match model {
ResolvedProductionModel::ConstantProductivity { productivity }
if *productivity > 0.0 =>
{
hb.max_turbined_m3s.min(hb.max_generation_mw / productivity)
}
_ => hb.max_turbined_m3s,
};
for blk in 0..layout.n_blks {
let col = layout.col_turbine_start + h_idx * layout.n_blks + blk;
col_lower[col] = hb.min_turbined_m3s;
col_upper[col] = turb_upper;
if is_fpha {
let block_hours = stage.blocks[blk].duration_hours;
objective[col] = hp.fpha_turbined_cost * block_hours;
}
}
}
for (h_idx, _hydro) in ctx.hydros.iter().enumerate() {
let hp = ctx.penalties.hydro_penalties(h_idx, stage_idx);
for blk in 0..layout.n_blks {
let col = layout.col_spillage_start + h_idx * layout.n_blks + blk;
col_upper[col] = f64::INFINITY;
let block_hours = stage.blocks[blk].duration_hours;
objective[col] = hp.spillage_cost * block_hours;
}
}
for (h_idx, hydro) in ctx.hydros.iter().enumerate() {
let hp = ctx.penalties.hydro_penalties(h_idx, stage_idx);
let max_div = hydro.diversion.as_ref().map_or(0.0, |d| d.max_flow_m3s);
for blk in 0..layout.n_blks {
let col = layout.col_diversion_start + h_idx * layout.n_blks + blk;
col_lower[col] = 0.0;
col_upper[col] = max_div;
if max_div > 0.0 {
let block_hours = stage.blocks[blk].duration_hours;
objective[col] = hp.diversion_cost * block_hours;
}
}
}
for (t_idx, thermal) in ctx.thermals.iter().enumerate() {
let tb = ctx.bounds.thermal_bounds(t_idx, stage_idx);
let marginal_cost_per_mwh = thermal
.cost_segments
.first()
.map_or(0.0, |seg| seg.cost_per_mwh);
for blk in 0..layout.n_blks {
let col = layout.col_thermal_start + t_idx * layout.n_blks + blk;
col_lower[col] = tb.min_generation_mw;
col_upper[col] = tb.max_generation_mw;
let block_hours = stage.blocks[blk].duration_hours;
objective[col] = marginal_cost_per_mwh * block_hours;
}
}
for (l_idx, line) in ctx.lines.iter().enumerate() {
let lb = ctx.bounds.line_bounds(l_idx, stage_idx);
let lp = ctx.penalties.line_penalties(l_idx, stage_idx);
for blk in 0..layout.n_blks {
let (df, rf) = ctx.resolved_exchange_factors.factors(l_idx, stage_idx, blk);
let col_fwd = layout.col_line_fwd_start + l_idx * layout.n_blks + blk;
let col_rev = layout.col_line_rev_start + l_idx * layout.n_blks + blk;
col_upper[col_fwd] = lb.direct_mw * df;
col_upper[col_rev] = lb.reverse_mw * rf;
let block_hours = stage.blocks[blk].duration_hours;
objective[col_fwd] = lp.exchange_cost * block_hours;
objective[col_rev] = lp.exchange_cost * block_hours;
let _ = line;
}
}
for (b_idx, bus) in ctx.buses.iter().enumerate() {
let bp = ctx.penalties.bus_penalties(b_idx, stage_idx);
for (seg_idx, segment) in bus.deficit_segments.iter().enumerate() {
for blk in 0..layout.n_blks {
let col_def = layout.col_deficit_start
+ b_idx * layout.max_deficit_segments * layout.n_blks
+ seg_idx * layout.n_blks
+ blk;
let block_hours = stage.blocks[blk].duration_hours;
col_upper[col_def] = segment.depth_mw.unwrap_or(f64::INFINITY);
objective[col_def] = segment.cost_per_mwh * block_hours;
}
}
for blk in 0..layout.n_blks {
let col_exc = layout.col_excess_start + b_idx * layout.n_blks + blk;
let block_hours = stage.blocks[blk].duration_hours;
col_upper[col_exc] = f64::INFINITY;
objective[col_exc] = bp.excess_cost * block_hours;
}
}
if ctx.has_penalty {
let total_stage_hours: f64 = stage.blocks.iter().map(|b| b.duration_hours).sum();
let penalty_cost = ctx.inflow_method.penalty_cost().unwrap_or(0.0);
let obj_coeff = penalty_cost * total_stage_hours;
for h_idx in 0..layout.n_h {
let col = layout.col_inflow_slack_start + h_idx;
objective[col] = obj_coeff;
}
}
for (local_idx, &h_idx) in layout.fpha_hydro_indices.iter().enumerate() {
let hb = ctx.bounds.hydro_bounds(h_idx, stage_idx);
for blk in 0..layout.n_blks {
let col = layout.col_generation_start + local_idx * layout.n_blks + blk;
col_lower[col] = 0.0;
col_upper[col] = hb.max_generation_mw;
}
}
let total_stage_hours: f64 = stage.blocks.iter().map(|b| b.duration_hours).sum();
for (local_idx, &h_idx) in layout.evap_hydro_indices.iter().enumerate() {
let col_q_ev = layout.col_evap_start + local_idx * 3;
let col_f_plus = layout.col_evap_start + local_idx * 3 + 1;
let col_f_minus = layout.col_evap_start + local_idx * 3 + 2;
col_lower[col_q_ev] = 0.0;
if let EvaporationModel::Linearized { coefficients, .. } =
ctx.evaporation_models.model(h_idx)
{
let coeff = &coefficients[stage_idx];
let hb = ctx.bounds.hydro_bounds(h_idx, stage_idx);
let q_ev_max = (coeff.k_evap0 + coeff.k_evap_v * hb.max_storage_hm3).max(0.0);
col_upper[col_q_ev] = q_ev_max * Q_EV_SAFETY_MARGIN;
} else {
col_upper[col_q_ev] = 0.0;
}
col_lower[col_f_plus] = 0.0;
col_upper[col_f_plus] = f64::INFINITY;
col_lower[col_f_minus] = 0.0;
col_upper[col_f_minus] = f64::INFINITY;
let hp = ctx.penalties.hydro_penalties(h_idx, stage_idx);
let obj = hp.evaporation_violation_cost * total_stage_hours;
objective[col_f_plus] = obj;
objective[col_f_minus] = obj * OVER_EVAPORATION_COST_MULTIPLIER;
}
for h_idx in 0..layout.n_h {
let col = layout.col_withdrawal_slack_start + h_idx;
let hb = ctx.bounds.hydro_bounds(h_idx, stage_idx);
if hb.water_withdrawal_m3s > 0.0 {
col_upper[col] = f64::INFINITY;
} else {
col_upper[col] = 0.0;
}
let hp = ctx.penalties.hydro_penalties(h_idx, stage_idx);
objective[col] = hp.water_withdrawal_violation_cost * total_stage_hours;
}
for (ncs_local, &ncs_sys_idx) in layout.active_ncs_indices.iter().enumerate() {
let ncs = &ctx.non_controllable_sources[ncs_sys_idx];
let avail_gen = ctx
.resolved_ncs_bounds
.available_generation(ncs_sys_idx, stage_idx);
for blk in 0..layout.n_blks {
let col = layout.col_ncs_start + ncs_local * layout.n_blks + blk;
let factor = ctx.resolved_ncs_factors.factor(ncs_sys_idx, stage_idx, blk);
col_upper[col] = avail_gen * factor;
let block_hours = stage.blocks[blk].duration_hours;
objective[col] = -ncs.curtailment_cost * block_hours;
}
}
for h_idx in 0..layout.n_h {
let col = layout.col_z_inflow_start + h_idx;
col_lower[col] = f64::NEG_INFINITY;
col_upper[col] = f64::INFINITY;
}
(col_lower, col_upper, objective)
}
fn fill_stage_rows(
ctx: &TemplateBuildCtx<'_>,
stage: &Stage,
stage_idx: usize,
layout: &StageLayout,
) -> (Vec<f64>, Vec<f64>) {
let mut row_lower = vec![0.0_f64; layout.num_rows];
let mut row_upper = vec![0.0_f64; layout.num_rows];
for h_idx in 0..layout.n_h {
let row = layout.row_water_balance_start + h_idx;
let base = if ctx.par_lp.n_stages() > 0 && ctx.par_lp.n_hydros() == layout.n_h {
ctx.par_lp.deterministic_base(stage_idx, h_idx)
} else {
0.0
};
let withdrawal = ctx
.bounds
.hydro_bounds(h_idx, stage_idx)
.water_withdrawal_m3s;
let rhs = layout.zeta * (base - withdrawal);
row_lower[row] = rhs;
row_upper[row] = rhs;
}
for (b_idx, bus) in ctx.buses.iter().enumerate() {
let mean_mw = ctx
.load_models
.iter()
.find(|lm| lm.bus_id == bus.id && lm.stage_id == stage.id)
.map_or(0.0, |lm| lm.mean_mw);
for blk in 0..layout.n_blks {
let factor = ctx.resolved_load_factors.factor(b_idx, stage_idx, blk);
let row = layout.row_load_balance_start + b_idx * layout.n_blks + blk;
let rhs = mean_mw * factor;
row_lower[row] = rhs;
row_upper[row] = rhs;
}
}
let n_blks = layout.n_blks;
for (local_idx, &h_idx) in layout.fpha_hydro_indices.iter().enumerate() {
if let ResolvedProductionModel::Fpha { planes, .. } =
ctx.production_models.model(h_idx, stage_idx)
{
let n_planes = planes.len();
debug_assert_eq!(
n_planes, layout.fpha_planes_per_hydro[local_idx],
"plane count mismatch for FPHA hydro {h_idx} at stage {stage_idx}"
);
for blk in 0..n_blks {
for (p_idx, plane) in planes.iter().enumerate() {
let row = layout.row_fpha_start
+ local_idx * n_blks * n_planes
+ blk * n_planes
+ p_idx;
row_lower[row] = f64::NEG_INFINITY;
row_upper[row] = plane.intercept;
}
}
}
}
for (local_idx, &h_idx) in layout.evap_hydro_indices.iter().enumerate() {
if let EvaporationModel::Linearized { coefficients, .. } =
ctx.evaporation_models.model(h_idx)
{
debug_assert!(
stage_idx < coefficients.len(),
"stage index {stage_idx} out of bounds for evaporation coefficients (len = {})",
coefficients.len()
);
let k_evap0 = coefficients[stage_idx].k_evap0;
let row = layout.row_evap_start + local_idx;
row_lower[row] = k_evap0;
row_upper[row] = k_evap0;
}
}
for h_idx in 0..layout.n_h {
let row = layout.row_z_inflow_start + h_idx;
let base = if ctx.par_lp.n_stages() > 0 && ctx.par_lp.n_hydros() == layout.n_h {
ctx.par_lp.deterministic_base(stage_idx, h_idx)
} else {
0.0
};
row_lower[row] = base;
row_upper[row] = base;
}
(row_lower, row_upper)
}
fn fill_state_and_water_entries(
ctx: &TemplateBuildCtx<'_>,
stage: &Stage,
stage_idx: usize,
layout: &StageLayout,
col_entries: &mut [Vec<(usize, f64)>],
) {
let idx = StageIndexer::new(ctx.n_hydros, ctx.max_par_order);
let n_h = layout.n_h;
let n_blks = layout.n_blks;
let lag_order = layout.lag_order;
let zeta = layout.zeta;
let row_water = layout.row_water_balance_start;
for h in 0..n_h {
let col = idx.storage_in.start + h;
col_entries[col].push((h, 1.0));
}
for lag in 0..lag_order {
for h in 0..n_h {
let col = idx.inflow_lags.start + lag * n_h + h;
let row = n_h + lag * n_h + h;
col_entries[col].push((row, 1.0));
}
}
for h_idx in 0..n_h {
let hydro = &ctx.hydros[h_idx];
let row = row_water + h_idx;
col_entries[h_idx].push((row, 1.0));
col_entries[idx.storage_in.start + h_idx].push((row, -1.0));
for blk in 0..n_blks {
let tau_h = stage.blocks[blk].duration_hours * M3S_TO_HM3;
let col_turbine = layout.col_turbine_start + h_idx * n_blks + blk;
col_entries[col_turbine].push((row, tau_h));
let col_spillage = layout.col_spillage_start + h_idx * n_blks + blk;
col_entries[col_spillage].push((row, tau_h));
let col_diversion = layout.col_diversion_start + h_idx * n_blks + blk;
col_entries[col_diversion].push((row, tau_h));
for &up_id in ctx.cascade.upstream(hydro.id) {
if let Some(&u_idx) = ctx.hydro_pos.get(&up_id) {
col_entries[layout.col_turbine_start + u_idx * n_blks + blk]
.push((row, -tau_h));
col_entries[layout.col_spillage_start + u_idx * n_blks + blk]
.push((row, -tau_h));
}
}
if let Some(sources) = ctx.diversion_upstream.get(&hydro.id) {
for &d_idx in sources {
let col_div = layout.col_diversion_start + d_idx * n_blks + blk;
col_entries[col_div].push((row, -tau_h));
}
}
}
if ctx.par_lp.n_stages() > 0 && ctx.par_lp.n_hydros() == n_h {
let psi = ctx.par_lp.psi_slice(stage_idx, h_idx);
for (lag, &psi_val) in psi.iter().enumerate() {
if psi_val != 0.0 && lag < lag_order {
let col = idx.inflow_lags.start + lag * n_h + h_idx;
col_entries[col].push((row, -zeta * psi_val));
}
}
}
}
if ctx.has_penalty {
for h_idx in 0..n_h {
let col = layout.col_inflow_slack_start + h_idx;
let row = row_water + h_idx;
col_entries[col].push((row, -zeta));
}
}
for (local_idx, &h_idx) in layout.evap_hydro_indices.iter().enumerate() {
let col_q_ev = layout.col_evap_start + local_idx * 3;
let row = row_water + h_idx;
col_entries[col_q_ev].push((row, zeta));
}
for h_idx in 0..n_h {
let col = layout.col_withdrawal_slack_start + h_idx;
let row = row_water + h_idx;
col_entries[col].push((row, -zeta));
}
}
fn fill_load_balance_entries(
ctx: &TemplateBuildCtx<'_>,
stage_idx: usize,
layout: &StageLayout,
col_entries: &mut [Vec<(usize, f64)>],
) {
let n_blks = layout.n_blks;
let row_load = layout.row_load_balance_start;
let mut fpha_local: Vec<Option<usize>> = vec![None; ctx.n_hydros];
for (local_idx, &h_idx) in layout.fpha_hydro_indices.iter().enumerate() {
fpha_local[h_idx] = Some(local_idx);
}
for (h_idx, hydro) in ctx.hydros.iter().enumerate() {
if let Some(local_idx) = fpha_local[h_idx] {
debug_assert!(
matches!(
ctx.production_models.model(h_idx, stage_idx),
ResolvedProductionModel::Fpha { .. }
),
"FPHA local-index table inconsistent with production model for hydro {h_idx}"
);
if let Some(&b_idx) = ctx.bus_pos.get(&hydro.bus_id) {
for blk in 0..n_blks {
let row = row_load + b_idx * n_blks + blk;
let col = layout.col_generation_start + local_idx * n_blks + blk;
col_entries[col].push((row, 1.0));
}
}
} else {
let rho = match &hydro.generation_model {
HydroGenerationModel::ConstantProductivity {
productivity_mw_per_m3s,
}
| HydroGenerationModel::LinearizedHead {
productivity_mw_per_m3s,
} => *productivity_mw_per_m3s,
HydroGenerationModel::Fpha => {
if let ResolvedProductionModel::ConstantProductivity { productivity } =
ctx.production_models.model(h_idx, stage_idx)
{
*productivity
} else {
debug_assert!(
false,
"Fpha entity model with non-Fpha resolved model and no local index for hydro {h_idx}"
);
0.0
}
}
};
if let Some(&b_idx) = ctx.bus_pos.get(&hydro.bus_id) {
for blk in 0..n_blks {
let row = row_load + b_idx * n_blks + blk;
let col = layout.col_turbine_start + h_idx * n_blks + blk;
col_entries[col].push((row, rho));
}
}
}
}
for (t_idx, thermal) in ctx.thermals.iter().enumerate() {
if let Some(&b_idx) = ctx.bus_pos.get(&thermal.bus_id) {
for blk in 0..n_blks {
let row = row_load + b_idx * n_blks + blk;
let col = layout.col_thermal_start + t_idx * n_blks + blk;
col_entries[col].push((row, 1.0));
}
}
}
for (l_idx, line) in ctx.lines.iter().enumerate() {
let src_idx = ctx.bus_pos.get(&line.source_bus_id).copied();
let tgt_idx = ctx.bus_pos.get(&line.target_bus_id).copied();
for blk in 0..n_blks {
let col_fwd = layout.col_line_fwd_start + l_idx * n_blks + blk;
let col_rev = layout.col_line_rev_start + l_idx * n_blks + blk;
if let Some(tgt) = tgt_idx {
let row = row_load + tgt * n_blks + blk;
col_entries[col_fwd].push((row, 1.0));
col_entries[col_rev].push((row, -1.0));
}
if let Some(src) = src_idx {
let row = row_load + src * n_blks + blk;
col_entries[col_fwd].push((row, -1.0));
col_entries[col_rev].push((row, 1.0));
}
}
}
for (b_idx, bus) in ctx.buses.iter().enumerate() {
for blk in 0..n_blks {
let row = row_load + b_idx * n_blks + blk;
for seg_idx in 0..bus.deficit_segments.len() {
let col_def = layout.col_deficit_start
+ b_idx * layout.max_deficit_segments * n_blks
+ seg_idx * n_blks
+ blk;
col_entries[col_def].push((row, 1.0));
}
let col_exc = layout.col_excess_start + b_idx * n_blks + blk;
col_entries[col_exc].push((row, -1.0));
}
}
}
fn fill_fpha_entries(
ctx: &TemplateBuildCtx<'_>,
stage_idx: usize,
layout: &StageLayout,
col_entries: &mut [Vec<(usize, f64)>],
) {
let idx = StageIndexer::new(ctx.n_hydros, ctx.max_par_order);
let n_blks = layout.n_blks;
for (local_idx, &h_idx) in layout.fpha_hydro_indices.iter().enumerate() {
let model = ctx.production_models.model(h_idx, stage_idx);
let planes = match model {
ResolvedProductionModel::Fpha { planes, .. } => planes,
ResolvedProductionModel::ConstantProductivity { .. } => {
debug_assert!(
false,
"fpha_hydro_indices contains hydro {h_idx} but model is ConstantProductivity"
);
continue;
}
};
let n_planes = planes.len();
for blk in 0..n_blks {
let col_v = h_idx; let col_v_in = idx.storage_in.start + h_idx; let col_q = layout.col_turbine_start + h_idx * n_blks + blk;
let col_s = layout.col_spillage_start + h_idx * n_blks + blk;
let col_g = layout.col_generation_start + local_idx * n_blks + blk;
for (p_idx, plane) in planes.iter().enumerate() {
let row =
layout.row_fpha_start + local_idx * n_blks * n_planes + blk * n_planes + p_idx;
col_entries[col_g].push((row, 1.0));
col_entries[col_v].push((row, -plane.gamma_v / 2.0));
col_entries[col_v_in].push((row, -plane.gamma_v / 2.0));
col_entries[col_q].push((row, -plane.gamma_q));
col_entries[col_s].push((row, -plane.gamma_s));
}
}
}
}
fn fill_evaporation_entries(
ctx: &TemplateBuildCtx<'_>,
stage_idx: usize,
layout: &StageLayout,
col_entries: &mut [Vec<(usize, f64)>],
) {
let idx = StageIndexer::new(ctx.n_hydros, ctx.max_par_order);
for (local_idx, &h_idx) in layout.evap_hydro_indices.iter().enumerate() {
let coeff = match ctx.evaporation_models.model(h_idx) {
EvaporationModel::Linearized { coefficients, .. } => {
debug_assert!(
stage_idx < coefficients.len(),
"evap_hydro_indices contains hydro {h_idx} but coefficients length {} \
is less than stage_idx {}",
coefficients.len(),
stage_idx
);
match coefficients.get(stage_idx) {
Some(c) => *c,
None => continue,
}
}
EvaporationModel::None => {
debug_assert!(
false,
"evap_hydro_indices contains hydro {h_idx} but model is None"
);
continue;
}
};
let col_q_ev = layout.col_evap_start + local_idx * 3;
let col_f_plus = layout.col_evap_start + local_idx * 3 + 1;
let col_f_minus = layout.col_evap_start + local_idx * 3 + 2;
let col_v = h_idx; let col_v_in = idx.storage_in.start + h_idx;
let row = layout.row_evap_start + local_idx;
col_entries[col_q_ev].push((row, 1.0));
col_entries[col_v].push((row, -coeff.k_evap_v / 2.0));
col_entries[col_v_in].push((row, -coeff.k_evap_v / 2.0));
col_entries[col_f_plus].push((row, 1.0));
col_entries[col_f_minus].push((row, -1.0));
}
}
#[allow(clippy::too_many_arguments)]
fn fill_generic_constraint_entries(
ctx: &TemplateBuildCtx<'_>,
stage: &Stage,
stage_idx: usize,
layout: &StageLayout,
col_entries: &mut [Vec<(usize, f64)>],
_col_lower: &mut [f64],
col_upper: &mut [f64],
objective: &mut [f64],
row_lower: &mut [f64],
row_upper: &mut [f64],
) {
if layout.n_generic_rows == 0 {
return;
}
let indexer = crate::indexer::StageIndexer::with_equipment_and_evaporation(
ctx.n_hydros,
ctx.max_par_order,
ctx.n_thermals,
ctx.n_lines,
ctx.n_buses,
layout.n_blks,
ctx.has_penalty,
layout.fpha_hydro_indices.clone(),
&layout.fpha_planes_per_hydro,
layout.evap_hydro_indices.clone(),
layout.max_deficit_segments,
);
for (entry_idx, entry) in layout.generic_constraint_rows.iter().enumerate() {
let row = layout.row_generic_start + entry_idx;
let constraint = &ctx.generic_constraints[entry.constraint_idx];
let block_hours = stage.blocks[entry.block_idx].duration_hours;
match entry.sense {
ConstraintSense::LessEqual => {
row_lower[row] = f64::NEG_INFINITY;
row_upper[row] = entry.bound;
}
ConstraintSense::GreaterEqual => {
row_lower[row] = entry.bound;
row_upper[row] = f64::INFINITY;
}
ConstraintSense::Equal => {
row_lower[row] = entry.bound;
row_upper[row] = entry.bound;
}
}
for term in &constraint.expression.terms {
let pairs = resolve_variable_ref(
&term.variable,
entry.block_idx,
layout.n_blks,
stage_idx,
&indexer,
ctx.production_models,
&ctx.hydro_pos,
&ctx.thermal_pos,
&ctx.bus_pos,
&ctx.line_pos,
);
for (col, multiplier) in pairs {
col_entries[col].push((row, term.coefficient * multiplier));
}
}
if let Some(plus_col) = entry.slack_plus_col {
let penalty = constraint.slack.penalty.unwrap_or(0.0);
let obj_coeff = penalty * block_hours;
col_upper[plus_col] = f64::INFINITY;
objective[plus_col] = obj_coeff;
match entry.sense {
ConstraintSense::LessEqual => {
col_entries[plus_col].push((row, -1.0));
}
ConstraintSense::GreaterEqual => {
col_entries[plus_col].push((row, 1.0));
}
ConstraintSense::Equal => {
col_entries[plus_col].push((row, 1.0));
}
}
if let Some(minus_col) = entry.slack_minus_col {
col_upper[minus_col] = f64::INFINITY;
objective[minus_col] = obj_coeff;
col_entries[minus_col].push((row, -1.0));
}
}
}
}
fn fill_ncs_load_balance_entries(
ctx: &TemplateBuildCtx<'_>,
layout: &StageLayout,
col_entries: &mut [Vec<(usize, f64)>],
) {
for (ncs_local, &ncs_sys_idx) in layout.active_ncs_indices.iter().enumerate() {
let ncs = &ctx.non_controllable_sources[ncs_sys_idx];
let Some(&bus_idx) = ctx.bus_pos.get(&ncs.bus_id) else {
continue;
};
for blk in 0..layout.n_blks {
let col = layout.col_ncs_start + ncs_local * layout.n_blks + blk;
let row = layout.row_load_balance_start + bus_idx * layout.n_blks + blk;
col_entries[col].push((row, 1.0));
}
}
}
fn fill_z_inflow_entries(
ctx: &TemplateBuildCtx<'_>,
stage_idx: usize,
layout: &StageLayout,
col_entries: &mut [Vec<(usize, f64)>],
) {
let n_h = layout.n_h;
let lag_order = layout.lag_order;
let idx = StageIndexer::new(n_h, lag_order);
for h_idx in 0..n_h {
let row = layout.row_z_inflow_start + h_idx;
let col_z = layout.col_z_inflow_start + h_idx;
col_entries[col_z].push((row, 1.0));
if ctx.par_lp.n_stages() > 0 && ctx.par_lp.n_hydros() == n_h {
let psi = ctx.par_lp.psi_slice(stage_idx, h_idx);
for (lag, &psi_val) in psi.iter().enumerate() {
if psi_val != 0.0 && lag < lag_order {
let col = idx.inflow_lags.start + lag * n_h + h_idx;
col_entries[col].push((row, -psi_val));
}
}
}
}
}
fn build_stage_matrix_entries(
ctx: &TemplateBuildCtx<'_>,
stage: &Stage,
stage_idx: usize,
layout: &StageLayout,
) -> Vec<Vec<(usize, f64)>> {
let mut col_entries: Vec<Vec<(usize, f64)>> = vec![Vec::new(); layout.num_cols];
fill_state_and_water_entries(ctx, stage, stage_idx, layout, &mut col_entries);
fill_load_balance_entries(ctx, stage_idx, layout, &mut col_entries);
fill_ncs_load_balance_entries(ctx, layout, &mut col_entries);
fill_fpha_entries(ctx, stage_idx, layout, &mut col_entries);
fill_evaporation_entries(ctx, stage_idx, layout, &mut col_entries);
fill_z_inflow_entries(ctx, stage_idx, layout, &mut col_entries);
col_entries
}
fn assemble_csc(col_entries: &[Vec<(usize, f64)>]) -> (Vec<i32>, Vec<i32>, Vec<f64>) {
let num_cols = col_entries.len();
let total_nz: usize = col_entries.iter().map(Vec::len).sum();
let mut col_starts = Vec::with_capacity(num_cols + 1);
let mut row_indices = Vec::with_capacity(total_nz);
let mut values = Vec::with_capacity(total_nz);
let mut offset: i32 = 0;
for entries in col_entries {
col_starts.push(offset);
for &(row, val) in entries {
#[allow(clippy::cast_possible_truncation, clippy::cast_possible_wrap)]
row_indices.push(row as i32);
values.push(val);
}
#[allow(clippy::cast_possible_truncation, clippy::cast_possible_wrap)]
{
offset += entries.len() as i32;
}
}
col_starts.push(offset);
(col_starts, row_indices, values)
}
#[allow(clippy::type_complexity)]
fn build_single_stage_template(
ctx: &TemplateBuildCtx<'_>,
stage: &Stage,
stage_idx: usize,
) -> (
StageTemplate,
usize,
usize,
Vec<GenericConstraintRowEntry>,
usize,
usize,
Vec<usize>,
) {
let layout = StageLayout::new(ctx, stage, stage_idx);
let stage_base_row = layout.row_water_balance_start;
let load_balance_row_start = layout.row_load_balance_start;
let (mut col_lower, mut col_upper, mut objective) =
fill_stage_columns(ctx, stage, stage_idx, &layout);
let (mut row_lower, mut row_upper) = fill_stage_rows(ctx, stage, stage_idx, &layout);
let mut col_entries = build_stage_matrix_entries(ctx, stage, stage_idx, &layout);
fill_generic_constraint_entries(
ctx,
stage,
stage_idx,
&layout,
&mut col_entries,
&mut col_lower,
&mut col_upper,
&mut objective,
&mut row_lower,
&mut row_upper,
);
let theta_col = StageIndexer::new(ctx.n_hydros, ctx.max_par_order).theta;
for (i, coeff) in objective.iter_mut().enumerate() {
if i != theta_col {
*coeff /= COST_SCALE_FACTOR;
}
}
for entries in &mut col_entries {
entries.sort_unstable_by_key(|&(row, _)| row);
}
let (col_starts, row_indices, values) = assemble_csc(&col_entries);
let n_transfer = ctx.n_hydros * ctx.max_par_order;
let total_nz = col_entries.iter().map(Vec::len).sum();
let gc_row_entries = layout.generic_constraint_rows;
let ncs_col_start = layout.col_ncs_start;
let n_ncs = layout.n_ncs;
let ncs_active = layout.active_ncs_indices;
let template = StageTemplate {
num_cols: layout.num_cols,
num_rows: layout.num_rows,
num_nz: total_nz,
col_starts,
row_indices,
values,
col_lower,
col_upper,
objective,
row_lower,
row_upper,
n_state: layout.n_state,
n_transfer,
n_dual_relevant: layout.n_dual_relevant,
n_hydro: layout.n_h,
max_par_order: layout.lag_order,
col_scale: Vec::new(),
row_scale: Vec::new(),
};
(
template,
stage_base_row,
load_balance_row_start,
gc_row_entries,
ncs_col_start,
n_ncs,
ncs_active,
)
}
#[must_use]
#[allow(clippy::cast_sign_loss)] pub(crate) fn compute_col_scale(num_cols: usize, col_starts: &[i32], values: &[f64]) -> Vec<f64> {
let mut scale = vec![1.0_f64; num_cols];
for j in 0..num_cols {
let start = col_starts[j] as usize;
let end = col_starts[j + 1] as usize;
if start == end {
continue;
}
let mut max_abs = 0.0_f64;
let mut min_abs = f64::INFINITY;
for &v in &values[start..end] {
let abs_val = v.abs();
if abs_val > 0.0 {
max_abs = max_abs.max(abs_val);
min_abs = min_abs.min(abs_val);
}
}
if max_abs > 0.0 && min_abs < f64::INFINITY {
let d = 1.0 / (max_abs * min_abs).sqrt();
scale[j] = d;
}
}
scale
}
pub(crate) fn apply_col_scale(template: &mut StageTemplate, col_scale: &[f64]) {
let num_cols = template.num_cols;
debug_assert_eq!(col_scale.len(), num_cols);
#[allow(clippy::needless_range_loop, clippy::cast_sign_loss)]
for j in 0..num_cols {
let start = template.col_starts[j] as usize;
let end = template.col_starts[j + 1] as usize;
let d = col_scale[j];
for v in &mut template.values[start..end] {
*v *= d;
}
}
for (obj, &d) in template.objective.iter_mut().zip(col_scale) {
*obj *= d;
}
for ((lo, hi), &d) in template
.col_lower
.iter_mut()
.zip(template.col_upper.iter_mut())
.zip(col_scale)
{
*lo /= d;
*hi /= d;
}
}
#[must_use]
#[allow(clippy::cast_sign_loss)] pub(crate) fn compute_row_scale(
num_rows: usize,
num_cols: usize,
col_starts: &[i32],
row_indices: &[i32],
values: &[f64],
) -> Vec<f64> {
let mut row_max = vec![0.0_f64; num_rows];
let mut row_min = vec![f64::INFINITY; num_rows];
#[allow(clippy::needless_range_loop)] for j in 0..num_cols {
let start = col_starts[j] as usize;
let end = col_starts[j + 1] as usize;
for k in start..end {
let row = row_indices[k] as usize;
let abs_val = values[k].abs();
if abs_val > 0.0 {
row_max[row] = row_max[row].max(abs_val);
row_min[row] = row_min[row].min(abs_val);
}
}
}
let mut scale = vec![1.0_f64; num_rows];
for (s, (&rmax, &rmin)) in scale.iter_mut().zip(row_max.iter().zip(row_min.iter())) {
if rmax > 0.0 && rmin < f64::INFINITY {
*s = 1.0 / (rmax * rmin).sqrt();
}
}
scale
}
pub(crate) fn apply_row_scale(template: &mut StageTemplate, row_scale: &[f64]) {
let num_rows = template.num_rows;
debug_assert_eq!(row_scale.len(), num_rows);
let num_cols = template.num_cols;
#[allow(clippy::needless_range_loop, clippy::cast_sign_loss)]
for j in 0..num_cols {
let start = template.col_starts[j] as usize;
let end = template.col_starts[j + 1] as usize;
for k in start..end {
let row = template.row_indices[k] as usize;
template.values[k] *= row_scale[row];
}
}
for ((lo, hi), &d) in template
.row_lower
.iter_mut()
.zip(template.row_upper.iter_mut())
.zip(row_scale)
{
*lo *= d;
*hi *= d;
}
}
fn compute_noise_scale(
study_stages: &[&Stage],
n_hydros: usize,
par_lp: &PrecomputedPar,
) -> (Vec<f64>, Vec<f64>, Vec<Vec<f64>>) {
let n = study_stages.len();
let mut noise_scale = vec![0.0_f64; n * n_hydros];
let mut zeta_per_stage = Vec::with_capacity(n);
let mut block_hours_per_stage = Vec::with_capacity(n);
for (s_idx, stage) in study_stages.iter().enumerate() {
let total_hours: f64 = stage.blocks.iter().map(|b| b.duration_hours).sum();
let zeta_s = total_hours * M3S_TO_HM3;
zeta_per_stage.push(zeta_s);
block_hours_per_stage.push(stage.blocks.iter().map(|b| b.duration_hours).collect());
for h_idx in 0..n_hydros {
let sigma = if par_lp.n_stages() > 0 && par_lp.n_hydros() == n_hydros {
par_lp.sigma(s_idx, h_idx)
} else {
0.0
};
noise_scale[s_idx * n_hydros + h_idx] = zeta_s * sigma;
}
}
(noise_scale, zeta_per_stage, block_hours_per_stage)
}
fn collect_load_bus_indices(system: &System, bus_pos: &HashMap<EntityId, usize>) -> Vec<usize> {
let mut ids: Vec<EntityId> = system
.load_models()
.iter()
.filter(|m| m.std_mw > 0.0)
.map(|m| m.bus_id)
.collect();
ids.sort_unstable_by_key(|id| id.0);
ids.dedup();
ids.iter()
.filter_map(|id| bus_pos.get(id).copied())
.collect()
}
#[allow(clippy::too_many_lines)]
pub fn build_stage_templates(
system: &System,
inflow_method: &InflowNonNegativityMethod,
par_lp: &PrecomputedPar,
normal_lp: &PrecomputedNormal,
production_models: &ProductionModelSet,
evaporation_models: &EvaporationModelSet,
) -> Result<StageTemplates, SddpError> {
let study_stages: Vec<_> = system.stages().iter().filter(|s| s.id >= 0).collect();
let hydros = system.hydros();
let n_hydros = hydros.len();
if study_stages.is_empty() {
return Ok(StageTemplates {
templates: Vec::new(),
base_rows: Vec::new(),
noise_scale: Vec::new(),
zeta_per_stage: Vec::new(),
block_hours_per_stage: Vec::new(),
n_hydros,
load_balance_row_starts: Vec::new(),
n_load_buses: 0,
load_bus_indices: Vec::new(),
generic_constraint_row_entries: Vec::new(),
ncs_col_starts: Vec::new(),
n_ncs_per_stage: Vec::new(),
active_ncs_indices: Vec::new(),
diversion_upstream: HashMap::new(),
hydro_productivities_per_stage: Vec::new(),
});
}
let buses = system.buses();
let hydro_pos: HashMap<EntityId, usize> =
hydros.iter().enumerate().map(|(i, h)| (h.id, i)).collect();
let thermal_pos: HashMap<EntityId, usize> = system
.thermals()
.iter()
.enumerate()
.map(|(i, t)| (t.id, i))
.collect();
let line_pos: HashMap<EntityId, usize> = system
.lines()
.iter()
.enumerate()
.map(|(i, l)| (l.id, i))
.collect();
let bus_pos: HashMap<EntityId, usize> =
buses.iter().enumerate().map(|(i, b)| (b.id, i)).collect();
let load_bus_indices = collect_load_bus_indices(system, &bus_pos);
let n_load_buses = load_bus_indices.len();
debug_assert!(
normal_lp.n_entities() == 0 || normal_lp.n_entities() == n_load_buses,
"PrecomputedNormal has {} entities but system has {} stochastic load buses",
normal_lp.n_entities(),
n_load_buses
);
let max_par_order: usize = system
.inflow_models()
.iter()
.filter(|m| m.stage_id >= 0)
.map(|m| m.ar_coefficients.len())
.max()
.unwrap_or(0);
let mut diversion_upstream: HashMap<EntityId, Vec<usize>> = HashMap::new();
for (h_idx, hydro) in hydros.iter().enumerate() {
if let Some(ref div) = hydro.diversion {
diversion_upstream
.entry(div.downstream_id)
.or_default()
.push(h_idx);
}
}
let diversion_upstream_output = diversion_upstream.clone();
let ctx = TemplateBuildCtx {
hydros,
thermals: system.thermals(),
lines: system.lines(),
buses,
load_models: system.load_models(),
cascade: system.cascade(),
bounds: system.bounds(),
penalties: system.penalties(),
hydro_pos,
thermal_pos,
line_pos,
bus_pos,
inflow_method,
par_lp,
production_models,
evaporation_models,
generic_constraints: system.generic_constraints(),
resolved_generic_bounds: system.resolved_generic_bounds(),
resolved_load_factors: system.resolved_load_factors(),
resolved_exchange_factors: system.resolved_exchange_factors(),
non_controllable_sources: system.non_controllable_sources(),
resolved_ncs_bounds: system.resolved_ncs_bounds(),
resolved_ncs_factors: system.resolved_ncs_factors(),
diversion_upstream,
n_hydros,
n_thermals: system.thermals().len(),
n_lines: system.lines().len(),
n_buses: buses.len(),
max_par_order,
has_penalty: n_hydros > 0 && inflow_method.has_slack_columns(),
};
let n_study = study_stages.len();
let mut templates = Vec::with_capacity(n_study);
let mut base_rows = Vec::with_capacity(n_study);
let mut load_balance_row_starts = Vec::with_capacity(n_study);
let mut generic_constraint_row_entries = Vec::with_capacity(n_study);
let mut ncs_col_starts = Vec::with_capacity(n_study);
let mut n_ncs_per_stage = Vec::with_capacity(n_study);
let mut active_ncs_indices_per_stage = Vec::with_capacity(n_study);
for (stage_idx, stage) in study_stages.iter().enumerate() {
let (
template,
stage_base_row,
load_balance_row_start,
gc_entries,
ncs_col_start,
ncs_count,
ncs_active,
) = build_single_stage_template(&ctx, stage, stage_idx);
templates.push(template);
base_rows.push(stage_base_row);
load_balance_row_starts.push(load_balance_row_start);
generic_constraint_row_entries.push(gc_entries);
ncs_col_starts.push(ncs_col_start);
n_ncs_per_stage.push(ncs_count);
active_ncs_indices_per_stage.push(ncs_active);
}
let (noise_scale, zeta_per_stage, block_hours_per_stage) =
compute_noise_scale(&study_stages, n_hydros, par_lp);
let hydro_productivities_per_stage: Vec<Vec<f64>> = (0..n_study)
.map(|s| {
(0..n_hydros)
.map(|h| match ctx.production_models.model(h, s) {
ResolvedProductionModel::ConstantProductivity { productivity } => *productivity,
ResolvedProductionModel::Fpha { .. } => 0.0,
})
.collect()
})
.collect();
Ok(StageTemplates {
templates,
base_rows,
noise_scale,
zeta_per_stage,
block_hours_per_stage,
n_hydros,
load_balance_row_starts,
n_load_buses,
load_bus_indices,
generic_constraint_row_entries,
ncs_col_starts,
n_ncs_per_stage,
active_ncs_indices: active_ncs_indices_per_stage,
diversion_upstream: diversion_upstream_output,
hydro_productivities_per_stage,
})
}
#[cfg(test)]
#[allow(
clippy::doc_markdown,
clippy::too_many_lines,
clippy::cast_sign_loss,
clippy::cast_possible_truncation
)]
mod tests {
use super::{
COST_SCALE_FACTOR, OVER_EVAPORATION_COST_MULTIPLIER, PatchBuffer, Q_EV_SAFETY_MARGIN,
ar_dynamics_row_offset,
};
use crate::indexer::StageIndexer;
fn idx(n: usize, l: usize) -> StageIndexer {
StageIndexer::new(n, l)
}
#[test]
fn new_3_2_sizes_to_15() {
let buf = PatchBuffer::new(3, 2, 0, 0);
assert_eq!(buf.indices.len(), 15);
assert_eq!(buf.lower.len(), 15);
assert_eq!(buf.upper.len(), 15);
}
#[test]
fn new_160_12_sizes_to_2400() {
let buf = PatchBuffer::new(160, 12, 0, 0);
assert_eq!(buf.indices.len(), 2400);
assert_eq!(buf.lower.len(), 2400);
assert_eq!(buf.upper.len(), 2400);
}
#[test]
fn new_zero_lags_sizes_to_3n() {
let buf = PatchBuffer::new(5, 0, 0, 0);
assert_eq!(buf.indices.len(), 15); }
#[test]
fn new_zero_hydros_sizes_to_zero() {
let buf = PatchBuffer::new(0, 0, 0, 0);
assert_eq!(buf.indices.len(), 0);
}
#[test]
fn forward_patch_count_without_z_inflow_fill() {
let buf = PatchBuffer::new(3, 2, 0, 0);
assert_eq!(buf.forward_patch_count(), 12);
}
#[test]
fn state_patch_count_is_n_times_one_plus_l() {
let buf = PatchBuffer::new(3, 2, 0, 0);
assert_eq!(buf.state_patch_count(), 9);
}
#[test]
fn state_patch_count_zero_lags() {
let buf = PatchBuffer::new(4, 0, 0, 0);
assert_eq!(buf.state_patch_count(), 4);
}
#[test]
fn ar_dynamics_row_offset_adds_base_plus_hydro() {
assert_eq!(ar_dynamics_row_offset(100, 0), 100);
assert_eq!(ar_dynamics_row_offset(100, 1), 101);
assert_eq!(ar_dynamics_row_offset(100, 2), 102);
}
#[test]
fn ar_dynamics_row_offset_zero_base() {
assert_eq!(ar_dynamics_row_offset(0, 7), 7);
}
#[test]
fn fill_forward_patches_category1_indices() {
let mut buf = PatchBuffer::new(3, 2, 0, 0);
let state = [10.0, 20.0, 30.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0];
let noise = [0.1, 0.2, 0.3];
buf.fill_forward_patches(&idx(3, 2), &state, &noise, 50, &[]);
assert_eq!(buf.indices[0], 0);
assert_eq!(buf.indices[1], 1);
assert_eq!(buf.indices[2], 2);
}
#[test]
fn fill_forward_patches_category2_indices() {
let mut buf = PatchBuffer::new(3, 2, 0, 0);
let state = [10.0, 20.0, 30.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0];
let noise = [0.1, 0.2, 0.3];
buf.fill_forward_patches(&idx(3, 2), &state, &noise, 50, &[]);
assert_eq!(buf.indices[3], 3); assert_eq!(buf.indices[4], 4); assert_eq!(buf.indices[5], 5); assert_eq!(buf.indices[6], 6); assert_eq!(buf.indices[7], 7); assert_eq!(buf.indices[8], 8); }
#[test]
fn fill_forward_patches_category3_indices() {
let mut buf = PatchBuffer::new(3, 2, 0, 0);
let state = [10.0, 20.0, 30.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0];
let noise = [0.1, 0.2, 0.3];
buf.fill_forward_patches(&idx(3, 2), &state, &noise, 50, &[]);
assert_eq!(buf.indices[9], 50); assert_eq!(buf.indices[10], 51); assert_eq!(buf.indices[11], 52); }
#[test]
fn fill_forward_patches_category1_values() {
let mut buf = PatchBuffer::new(3, 2, 0, 0);
let state = [10.0, 20.0, 30.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0];
let noise = [0.1, 0.2, 0.3];
buf.fill_forward_patches(&idx(3, 2), &state, &noise, 50, &[]);
assert_eq!(buf.lower[0], 10.0);
assert_eq!(buf.upper[0], 10.0);
assert_eq!(buf.lower[1], 20.0);
assert_eq!(buf.upper[1], 20.0);
assert_eq!(buf.lower[2], 30.0);
assert_eq!(buf.upper[2], 30.0);
}
#[test]
fn fill_forward_patches_category2_values() {
let mut buf = PatchBuffer::new(3, 2, 0, 0);
let state = [10.0, 20.0, 30.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0];
let noise = [0.1, 0.2, 0.3];
buf.fill_forward_patches(&idx(3, 2), &state, &noise, 50, &[]);
assert_eq!(buf.lower[3], 1.0);
assert_eq!(buf.upper[3], 1.0);
assert_eq!(buf.lower[4], 2.0);
assert_eq!(buf.upper[4], 2.0);
assert_eq!(buf.lower[5], 3.0);
assert_eq!(buf.upper[5], 3.0);
assert_eq!(buf.lower[6], 4.0);
assert_eq!(buf.upper[6], 4.0);
assert_eq!(buf.lower[7], 5.0);
assert_eq!(buf.upper[7], 5.0);
assert_eq!(buf.lower[8], 6.0);
assert_eq!(buf.upper[8], 6.0);
}
#[test]
fn fill_forward_patches_category3_values() {
let mut buf = PatchBuffer::new(3, 2, 0, 0);
let state = [10.0, 20.0, 30.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0];
let noise = [0.1, 0.2, 0.3];
buf.fill_forward_patches(&idx(3, 2), &state, &noise, 50, &[]);
assert_eq!(buf.lower[9], 0.1);
assert_eq!(buf.upper[9], 0.1);
assert_eq!(buf.lower[10], 0.2);
assert_eq!(buf.upper[10], 0.2);
assert_eq!(buf.lower[11], 0.3);
assert_eq!(buf.upper[11], 0.3);
}
#[test]
fn fill_forward_patches_all_equality_constraints() {
let mut buf = PatchBuffer::new(3, 2, 0, 0);
let state = [10.0, 20.0, 30.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0];
let noise = [0.1, 0.2, 0.3];
buf.fill_forward_patches(&idx(3, 2), &state, &noise, 50, &[]);
for i in 0..buf.forward_patch_count() {
assert_eq!(
buf.lower[i],
buf.upper[i],
"patch {i}: lower {lo} != upper {up}",
lo = buf.lower[i],
up = buf.upper[i],
);
}
}
#[test]
fn fill_state_patches_count_is_n_times_one_plus_l() {
let mut buf = PatchBuffer::new(3, 2, 0, 0);
let state = [10.0, 20.0, 30.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0];
buf.fill_state_patches(&idx(3, 2), &state, &[]);
assert_eq!(buf.state_patch_count(), 9);
}
#[test]
fn fill_state_patches_category1_correct() {
let mut buf = PatchBuffer::new(3, 2, 0, 0);
let state = [10.0, 20.0, 30.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0];
buf.fill_state_patches(&idx(3, 2), &state, &[]);
assert_eq!(buf.indices[0], 0);
assert_eq!(buf.lower[0], 10.0);
assert_eq!(buf.upper[0], 10.0);
assert_eq!(buf.indices[1], 1);
assert_eq!(buf.lower[1], 20.0);
assert_eq!(buf.upper[1], 20.0);
assert_eq!(buf.indices[2], 2);
assert_eq!(buf.lower[2], 30.0);
assert_eq!(buf.upper[2], 30.0);
}
#[test]
fn fill_state_patches_category2_correct() {
let mut buf = PatchBuffer::new(3, 2, 0, 0);
let state = [10.0, 20.0, 30.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0];
buf.fill_state_patches(&idx(3, 2), &state, &[]);
assert_eq!(buf.indices[3], 3);
assert_eq!(buf.lower[3], 1.0);
assert_eq!(buf.upper[3], 1.0);
assert_eq!(buf.indices[8], 8);
assert_eq!(buf.lower[8], 6.0);
assert_eq!(buf.upper[8], 6.0);
}
#[test]
fn fill_state_patches_equality_constraints_in_active_range() {
let mut buf = PatchBuffer::new(3, 2, 0, 0);
let state = [10.0, 20.0, 30.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0];
buf.fill_state_patches(&idx(3, 2), &state, &[]);
let active = buf.state_patch_count();
for i in 0..active {
assert_eq!(
buf.lower[i],
buf.upper[i],
"state patch {i}: lower {lo} != upper {up}",
lo = buf.lower[i],
up = buf.upper[i],
);
}
}
#[test]
fn forward_patches_zero_lags_only_storage_and_noise() {
let n = 2;
let mut buf = PatchBuffer::new(n, 0, 0, 0);
let state = [5.0, 7.0]; let noise = [0.5, 0.6];
buf.fill_forward_patches(&idx(n, 0), &state, &noise, 10, &[]);
assert_eq!(buf.forward_patch_count(), 4);
assert_eq!(buf.indices[0], 0);
assert_eq!(buf.lower[0], 5.0);
assert_eq!(buf.indices[1], 1);
assert_eq!(buf.lower[1], 7.0);
assert_eq!(buf.indices[2], 10); assert_eq!(buf.lower[2], 0.5);
assert_eq!(buf.indices[3], 11); assert_eq!(buf.lower[3], 0.6);
}
#[test]
fn state_patches_zero_lags_only_storage() {
let n = 3;
let mut buf = PatchBuffer::new(n, 0, 0, 0);
let state = [1.0, 2.0, 3.0]; buf.fill_state_patches(&idx(n, 0), &state, &[]);
assert_eq!(buf.state_patch_count(), 3);
assert_eq!(buf.indices[0], 0);
assert_eq!(buf.lower[0], 1.0);
assert_eq!(buf.upper[0], 1.0);
assert_eq!(buf.indices[1], 1);
assert_eq!(buf.lower[1], 2.0);
assert_eq!(buf.indices[2], 2);
assert_eq!(buf.lower[2], 3.0);
}
#[test]
fn production_scale_forward_patch_count() {
let buf = PatchBuffer::new(160, 12, 0, 0);
assert_eq!(buf.forward_patch_count(), 2240);
assert_eq!(buf.indices.len(), 2400);
}
#[test]
#[allow(clippy::cast_precision_loss)] fn production_scale_fill_forward_patches_smoke() {
let n = 160;
let l = 12;
let mut buf = PatchBuffer::new(n, l, 0, 0);
let n_state = n * (1 + l);
let state: Vec<f64> = (0..n_state).map(|i| i as f64).collect();
let noise: Vec<f64> = (0..n).map(|h| h as f64 * 0.01).collect();
buf.fill_forward_patches(&StageIndexer::new(n, l), &state, &noise, 500, &[]);
assert_eq!(buf.indices[0], 0);
assert_eq!(buf.lower[0], 0.0);
assert_eq!(buf.indices[2080], 500); assert_eq!(buf.lower[2080], 0.0); assert_eq!(buf.indices[2239], 659); assert_eq!(buf.lower[2239], 159.0 * 0.01);
for i in 0..buf.forward_patch_count() {
assert_eq!(buf.lower[i], buf.upper[i], "patch {i} not equality");
}
}
#[test]
fn clone_and_debug() {
let buf = PatchBuffer::new(3, 2, 0, 0);
let cloned = buf.clone();
assert_eq!(cloned.indices.len(), buf.indices.len());
let s = format!("{buf:?}");
assert!(s.contains("PatchBuffer"));
}
#[test]
fn new_with_load_allocates_correct_capacity() {
let buf = PatchBuffer::new(2, 1, 1, 3);
assert_eq!(buf.indices.len(), 11);
assert_eq!(buf.lower.len(), 11);
assert_eq!(buf.upper.len(), 11);
}
#[test]
fn fill_load_patches_correct_indices() {
let mut buf = PatchBuffer::new(0, 0, 2, 2);
let load_rhs = [300.0_f64, 280.0, 500.0, 450.0];
let bus_positions = [0_usize, 1];
buf.fill_load_patches(100, 2, &load_rhs, &bus_positions, &[]);
assert_eq!(buf.indices[0], 100); assert_eq!(buf.indices[1], 101); assert_eq!(buf.indices[2], 102); assert_eq!(buf.indices[3], 103); }
#[test]
fn fill_load_patches_correct_values() {
let mut buf = PatchBuffer::new(0, 0, 2, 2);
let load_rhs = [300.0_f64, 280.0, 500.0, 450.0];
let bus_positions = [0_usize, 1];
buf.fill_load_patches(100, 2, &load_rhs, &bus_positions, &[]);
assert_eq!(buf.lower[0], 300.0);
assert_eq!(buf.upper[0], 300.0);
assert_eq!(buf.lower[1], 280.0);
assert_eq!(buf.upper[1], 280.0);
assert_eq!(buf.lower[2], 500.0);
assert_eq!(buf.upper[2], 500.0);
assert_eq!(buf.lower[3], 450.0);
assert_eq!(buf.upper[3], 450.0);
}
#[test]
fn fill_load_patches_equality_constraints() {
let mut buf = PatchBuffer::new(3, 2, 2, 3);
let state = [10.0, 20.0, 30.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0];
let noise = [0.1, 0.2, 0.3];
buf.fill_forward_patches(&idx(3, 2), &state, &noise, 50, &[]);
let load_rhs = [100.0_f64, 90.0, 80.0, 200.0, 190.0, 180.0];
let bus_positions = [0_usize, 1];
buf.fill_load_patches(20, 3, &load_rhs, &bus_positions, &[]);
let count = buf.forward_patch_count();
for i in 0..count {
assert_eq!(
buf.lower[i],
buf.upper[i],
"patch {i}: lower {lo} != upper {up}",
lo = buf.lower[i],
up = buf.upper[i],
);
}
}
#[test]
fn forward_patch_count_includes_load() {
let mut buf = PatchBuffer::new(3, 2, 2, 3);
let state = [10.0, 20.0, 30.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0];
let noise = [0.1, 0.2, 0.3];
buf.fill_forward_patches(&idx(3, 2), &state, &noise, 50, &[]);
let load_rhs = [100.0_f64, 90.0, 80.0, 200.0, 190.0, 180.0];
let bus_positions = [0_usize, 1];
buf.fill_load_patches(20, 3, &load_rhs, &bus_positions, &[]);
assert_eq!(buf.forward_patch_count(), 18); }
#[test]
fn state_patch_count_excludes_load() {
let mut buf = PatchBuffer::new(3, 2, 2, 3);
let state = [10.0, 20.0, 30.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0];
buf.fill_state_patches(&idx(3, 2), &state, &[]);
let load_rhs = [100.0_f64, 90.0, 80.0, 200.0, 190.0, 180.0];
let bus_positions = [0_usize, 1];
buf.fill_load_patches(20, 3, &load_rhs, &bus_positions, &[]);
assert_eq!(buf.state_patch_count(), 9);
}
#[test]
fn zero_load_buses_no_category4() {
let mut buf = PatchBuffer::new(3, 2, 0, 0);
let state = [10.0, 20.0, 30.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0];
let noise = [0.1, 0.2, 0.3];
buf.fill_forward_patches(&idx(3, 2), &state, &noise, 50, &[]);
assert_eq!(buf.forward_patch_count(), 12); }
use super::build_stage_templates;
use crate::hydro_models::{
FphaPlane, PrepareHydroModelsResult, ProductionModelSet, ResolvedProductionModel,
};
use crate::inflow_method::InflowNonNegativityMethod;
use cobre_core::{
Bus, BusStagePenalties, ContractStageBounds, DeficitSegment, EntityId, HydroStageBounds,
HydroStagePenalties, LineStageBounds, LineStagePenalties, NcsStagePenalties,
PumpingStageBounds, ResolvedBounds, ResolvedPenalties, SystemBuilder, ThermalStageBounds,
};
use cobre_stochastic::normal::precompute::PrecomputedNormal;
use cobre_stochastic::par::precompute::PrecomputedPar;
use crate::hydro_models::{EvaporationModel, EvaporationModelSet};
fn default_production(system: &cobre_core::System) -> ProductionModelSet {
PrepareHydroModelsResult::default_from_system(system).production
}
fn default_evaporation(system: &cobre_core::System) -> EvaporationModelSet {
PrepareHydroModelsResult::default_from_system(system).evaporation
}
fn no_penalty_config() -> InflowNonNegativityMethod {
InflowNonNegativityMethod::None
}
fn penalty_config(cost: f64) -> InflowNonNegativityMethod {
InflowNonNegativityMethod::Penalty { cost }
}
fn default_hydro_bounds() -> HydroStageBounds {
HydroStageBounds {
min_storage_hm3: 0.0,
max_storage_hm3: 200.0,
min_turbined_m3s: 0.0,
max_turbined_m3s: 100.0,
min_outflow_m3s: 0.0,
max_outflow_m3s: None,
min_generation_mw: 0.0,
max_generation_mw: 250.0,
max_diversion_m3s: None,
filling_inflow_m3s: 0.0,
water_withdrawal_m3s: 0.0,
}
}
fn default_hydro_penalties() -> HydroStagePenalties {
HydroStagePenalties {
spillage_cost: 0.01,
diversion_cost: 0.0,
fpha_turbined_cost: 0.0,
storage_violation_below_cost: 0.0,
filling_target_violation_cost: 0.0,
turbined_violation_below_cost: 0.0,
outflow_violation_below_cost: 0.0,
outflow_violation_above_cost: 0.0,
generation_violation_below_cost: 0.0,
evaporation_violation_cost: 0.0,
water_withdrawal_violation_cost: 0.0,
}
}
#[allow(clippy::cast_possible_truncation, clippy::cast_possible_wrap)]
fn one_bus_system(n_stages: usize) -> cobre_core::System {
use chrono::NaiveDate;
use cobre_core::scenario::LoadModel;
use cobre_core::temporal::{
Block, BlockMode, NoiseMethod, ScenarioSourceConfig, Stage, StageRiskConfig,
StageStateConfig,
};
let bus = Bus {
id: EntityId(1),
name: "B1".to_string(),
deficit_segments: vec![DeficitSegment {
depth_mw: None,
cost_per_mwh: 500.0,
}],
excess_cost: 0.0,
};
let stages: Vec<Stage> = (0..n_stages)
.map(|i| Stage {
index: i,
id: i as i32,
start_date: NaiveDate::from_ymd_opt(2024, 1, 1).unwrap(),
end_date: NaiveDate::from_ymd_opt(2024, 2, 1).unwrap(),
season_id: None,
blocks: vec![Block {
index: 0,
name: "S".to_string(),
duration_hours: 744.0,
}],
block_mode: BlockMode::Parallel,
state_config: StageStateConfig {
storage: false,
inflow_lags: false,
},
risk_config: StageRiskConfig::Expectation,
scenario_config: ScenarioSourceConfig {
branching_factor: 1,
noise_method: NoiseMethod::Saa,
},
})
.collect();
let load_models: Vec<LoadModel> = (0..n_stages)
.map(|i| LoadModel {
bus_id: EntityId(1),
stage_id: i as i32,
mean_mw: 100.0,
std_mw: 0.0,
})
.collect();
let n_st = n_stages.max(1);
let bounds = ResolvedBounds::new(
0,
0,
0,
0,
0,
n_st,
default_hydro_bounds(),
ThermalStageBounds {
min_generation_mw: 0.0,
max_generation_mw: 0.0,
},
LineStageBounds {
direct_mw: 0.0,
reverse_mw: 0.0,
},
PumpingStageBounds {
min_flow_m3s: 0.0,
max_flow_m3s: 0.0,
},
ContractStageBounds {
min_mw: 0.0,
max_mw: 0.0,
price_per_mwh: 0.0,
},
);
let penalties = ResolvedPenalties::new(
0,
1,
0,
0,
n_st,
default_hydro_penalties(),
BusStagePenalties { excess_cost: 0.0 },
LineStagePenalties { exchange_cost: 0.0 },
NcsStagePenalties {
curtailment_cost: 0.0,
},
);
SystemBuilder::new()
.buses(vec![bus])
.stages(stages)
.load_models(load_models)
.bounds(bounds)
.penalties(penalties)
.build()
.expect("one_bus_system: valid")
}
#[allow(
clippy::cast_possible_truncation,
clippy::cast_possible_wrap,
clippy::too_many_lines
)]
fn one_hydro_system(n_stages: usize, lag_order: usize) -> cobre_core::System {
use chrono::NaiveDate;
use cobre_core::entities::hydro::{Hydro, HydroGenerationModel, HydroPenalties};
use cobre_core::scenario::{InflowModel, LoadModel};
use cobre_core::temporal::{
Block, BlockMode, NoiseMethod, ScenarioSourceConfig, Stage, StageRiskConfig,
StageStateConfig,
};
let bus = Bus {
id: EntityId(1),
name: "B1".to_string(),
deficit_segments: vec![DeficitSegment {
depth_mw: None,
cost_per_mwh: 500.0,
}],
excess_cost: 0.0,
};
let hydro = Hydro {
id: EntityId(2),
name: "H1".to_string(),
bus_id: EntityId(1),
downstream_id: None,
entry_stage_id: None,
exit_stage_id: None,
min_storage_hm3: 0.0,
max_storage_hm3: 200.0,
min_outflow_m3s: 0.0,
max_outflow_m3s: None,
generation_model: HydroGenerationModel::ConstantProductivity {
productivity_mw_per_m3s: 2.5,
},
min_turbined_m3s: 0.0,
max_turbined_m3s: 100.0,
min_generation_mw: 0.0,
max_generation_mw: 250.0,
tailrace: None,
hydraulic_losses: None,
efficiency: None,
evaporation_coefficients_mm: None,
evaporation_reference_volumes_hm3: None,
diversion: None,
filling: None,
penalties: HydroPenalties {
spillage_cost: 0.01,
diversion_cost: 0.0,
fpha_turbined_cost: 0.0,
storage_violation_below_cost: 0.0,
filling_target_violation_cost: 0.0,
turbined_violation_below_cost: 0.0,
outflow_violation_below_cost: 0.0,
outflow_violation_above_cost: 0.0,
generation_violation_below_cost: 0.0,
evaporation_violation_cost: 0.0,
water_withdrawal_violation_cost: 0.0,
},
};
let stages: Vec<Stage> = (0..n_stages)
.map(|i| Stage {
index: i,
id: i as i32,
start_date: NaiveDate::from_ymd_opt(2024, 1, 1).unwrap(),
end_date: NaiveDate::from_ymd_opt(2024, 2, 1).unwrap(),
season_id: None,
blocks: vec![Block {
index: 0,
name: "S".to_string(),
duration_hours: 744.0,
}],
block_mode: BlockMode::Parallel,
state_config: StageStateConfig {
storage: true,
inflow_lags: lag_order > 0,
},
risk_config: StageRiskConfig::Expectation,
scenario_config: ScenarioSourceConfig {
branching_factor: 1,
noise_method: NoiseMethod::Saa,
},
})
.collect();
let ar_coefficients: Vec<f64> = (0..lag_order).map(|_| 0.5).collect();
let inflow_models: Vec<InflowModel> = (0..n_stages)
.map(|i| InflowModel {
hydro_id: EntityId(2),
stage_id: i as i32,
mean_m3s: 80.0,
std_m3s: 20.0,
ar_coefficients: ar_coefficients.clone(),
residual_std_ratio: 1.0,
})
.collect();
let load_models: Vec<LoadModel> = (0..n_stages)
.map(|i| LoadModel {
bus_id: EntityId(1),
stage_id: i as i32,
mean_mw: 100.0,
std_mw: 0.0,
})
.collect();
let n_st = n_stages.max(1);
let bounds = ResolvedBounds::new(
1,
0,
0,
0,
0,
n_st,
default_hydro_bounds(),
ThermalStageBounds {
min_generation_mw: 0.0,
max_generation_mw: 0.0,
},
LineStageBounds {
direct_mw: 0.0,
reverse_mw: 0.0,
},
PumpingStageBounds {
min_flow_m3s: 0.0,
max_flow_m3s: 0.0,
},
ContractStageBounds {
min_mw: 0.0,
max_mw: 0.0,
price_per_mwh: 0.0,
},
);
let penalties = ResolvedPenalties::new(
1,
1,
0,
0,
n_st,
default_hydro_penalties(),
BusStagePenalties { excess_cost: 0.0 },
LineStagePenalties { exchange_cost: 0.0 },
NcsStagePenalties {
curtailment_cost: 0.0,
},
);
SystemBuilder::new()
.buses(vec![bus])
.hydros(vec![hydro])
.stages(stages)
.inflow_models(inflow_models)
.load_models(load_models)
.bounds(bounds)
.penalties(penalties)
.build()
.expect("one_hydro_system: valid")
}
#[test]
fn empty_stages_returns_empty() {
let system = one_bus_system(0);
let result = build_stage_templates(
&system,
&no_penalty_config(),
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&default_production(&system),
&default_evaporation(&system),
)
.expect("constant productivity ok");
assert!(result.templates.is_empty());
assert!(result.base_rows.is_empty());
}
#[test]
fn one_stage_one_template() {
let system = one_bus_system(1);
let result = build_stage_templates(
&system,
&no_penalty_config(),
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&default_production(&system),
&default_evaporation(&system),
)
.expect("constant productivity ok");
assert_eq!(result.templates.len(), 1);
assert_eq!(result.base_rows.len(), 1);
}
#[test]
fn num_cols_formula_no_hydro_no_thermal_no_line() {
let system = one_bus_system(1);
let result = build_stage_templates(
&system,
&no_penalty_config(),
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&default_production(&system),
&default_evaporation(&system),
)
.expect("constant productivity ok");
let t = &result.templates[0];
assert_eq!(t.num_cols, 3, "num_cols mismatch for no-entity system");
}
#[test]
fn num_cols_formula_one_hydro_lag_zero() {
let system = one_hydro_system(1, 0);
let result = build_stage_templates(
&system,
&no_penalty_config(),
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&default_production(&system),
&default_evaporation(&system),
)
.expect("constant productivity ok");
let t = &result.templates[0];
assert_eq!(t.num_cols, 10, "num_cols mismatch for N=1 L=0");
}
#[test]
fn num_cols_formula_one_hydro_lag_two() {
let system = one_hydro_system(1, 2);
let result = build_stage_templates(
&system,
&no_penalty_config(),
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&default_production(&system),
&default_evaporation(&system),
)
.expect("constant productivity ok");
let t = &result.templates[0];
assert_eq!(t.num_cols, 12, "num_cols mismatch for N=1 L=2");
}
#[test]
fn num_rows_formula_no_hydro() {
let system = one_bus_system(1);
let result = build_stage_templates(
&system,
&no_penalty_config(),
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&default_production(&system),
&default_evaporation(&system),
)
.expect("constant productivity ok");
let t = &result.templates[0];
assert_eq!(t.num_rows, 1, "num_rows mismatch for no-hydro system");
}
#[test]
fn num_rows_formula_one_hydro_lag_zero() {
let system = one_hydro_system(1, 0);
let result = build_stage_templates(
&system,
&no_penalty_config(),
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&default_production(&system),
&default_evaporation(&system),
)
.expect("constant productivity ok");
let t = &result.templates[0];
assert_eq!(t.num_rows, 4, "num_rows mismatch for N=1 L=0");
}
#[test]
fn num_rows_formula_one_hydro_lag_two() {
let system = one_hydro_system(1, 2);
let result = build_stage_templates(
&system,
&no_penalty_config(),
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&default_production(&system),
&default_evaporation(&system),
)
.expect("constant productivity ok");
let t = &result.templates[0];
assert_eq!(t.num_rows, 6, "num_rows mismatch for N=1 L=2");
}
#[test]
fn n_state_matches_indexer() {
let system = one_hydro_system(1, 2);
let result = build_stage_templates(
&system,
&no_penalty_config(),
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&default_production(&system),
&default_evaporation(&system),
)
.expect("constant productivity ok");
let t = &result.templates[0];
let expected = StageIndexer::new(1, 2).n_state;
assert_eq!(t.n_state, expected, "n_state must match StageIndexer");
}
#[test]
fn n_transfer_is_n_times_lag_order() {
let system = one_hydro_system(1, 2);
let result = build_stage_templates(
&system,
&no_penalty_config(),
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&default_production(&system),
&default_evaporation(&system),
)
.expect("constant productivity ok");
let t = &result.templates[0];
assert_eq!(t.n_transfer, 2, "n_transfer = N*L");
}
#[test]
fn n_dual_relevant_equals_n_state_for_constant_productivity() {
let system = one_hydro_system(1, 2);
let result = build_stage_templates(
&system,
&no_penalty_config(),
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&default_production(&system),
&default_evaporation(&system),
)
.expect("constant productivity ok");
let t = &result.templates[0];
assert_eq!(
t.n_dual_relevant, t.n_state,
"n_dual_relevant must equal n_state for constant-productivity hydros"
);
}
#[test]
fn base_row_is_n_dual_relevant_plus_n_hydros() {
let system = one_hydro_system(2, 2);
let result = build_stage_templates(
&system,
&no_penalty_config(),
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&default_production(&system),
&default_evaporation(&system),
)
.expect("constant productivity ok");
for (s, (&br, t)) in result.base_rows.iter().zip(&result.templates).enumerate() {
assert_eq!(
br,
t.n_dual_relevant + t.n_hydro,
"base_rows[{s}] must equal n_dual_relevant + n_hydro"
);
}
}
#[test]
fn csc_col_starts_monotone_nondecreasing() {
let system = one_hydro_system(1, 1);
let result = build_stage_templates(
&system,
&no_penalty_config(),
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&default_production(&system),
&default_evaporation(&system),
)
.expect("constant productivity ok");
let t = &result.templates[0];
for w in t.col_starts.windows(2) {
assert!(w[0] <= w[1], "col_starts not monotone: {} > {}", w[0], w[1]);
}
assert_eq!(t.col_starts.len(), t.num_cols + 1);
}
#[test]
#[allow(clippy::cast_sign_loss)]
fn csc_row_indices_in_range() {
let system = one_hydro_system(1, 1);
let result = build_stage_templates(
&system,
&no_penalty_config(),
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&default_production(&system),
&default_evaporation(&system),
)
.expect("constant productivity ok");
let t = &result.templates[0];
for &r in &t.row_indices {
assert!(
r >= 0 && (r as usize) < t.num_rows,
"row index {r} out of range [0, {})",
t.num_rows
);
}
}
#[test]
#[allow(clippy::cast_sign_loss)]
fn csc_nz_count_matches_col_starts() {
let system = one_hydro_system(1, 1);
let result = build_stage_templates(
&system,
&no_penalty_config(),
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&default_production(&system),
&default_evaporation(&system),
)
.expect("constant productivity ok");
let t = &result.templates[0];
assert_eq!(
t.num_nz,
*t.col_starts.last().unwrap() as usize,
"num_nz must equal col_starts[num_cols]"
);
assert_eq!(
t.row_indices.len(),
t.num_nz,
"row_indices.len() must equal num_nz"
);
assert_eq!(t.values.len(), t.num_nz, "values.len() must equal num_nz");
}
#[test]
fn theta_column_has_unit_objective() {
let lag_order = 2;
let system = one_hydro_system(1, lag_order);
let result = build_stage_templates(
&system,
&no_penalty_config(),
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&default_production(&system),
&default_evaporation(&system),
)
.expect("constant productivity ok");
let t = &result.templates[0];
let theta_col = StageIndexer::new(1, lag_order).theta;
assert_eq!(
t.objective[theta_col], 1.0,
"theta column objective must be 1.0 (theta is not scaled by COST_SCALE_FACTOR)"
);
}
#[test]
fn spillage_objective_nonzero_for_nonzero_penalty() {
let system = one_hydro_system(1, 0);
let result = build_stage_templates(
&system,
&no_penalty_config(),
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&default_production(&system),
&default_evaporation(&system),
)
.expect("constant productivity ok");
let t = &result.templates[0];
let spill_col = 5;
assert!(
t.objective[spill_col] > 0.0,
"spillage objective must be > 0 when spillage_cost > 0"
);
}
#[allow(
clippy::cast_possible_truncation,
clippy::cast_possible_wrap,
clippy::too_many_lines
)]
fn fpha_system_with_turbined_cost(
n_planes: usize,
fpha_turbined_cost: f64,
block_durations_hours: &[f64],
) -> (cobre_core::System, ProductionModelSet) {
use chrono::NaiveDate;
use cobre_core::entities::hydro::{Hydro, HydroGenerationModel, HydroPenalties};
use cobre_core::scenario::{InflowModel, LoadModel};
use cobre_core::temporal::{
Block, BlockMode, NoiseMethod, ScenarioSourceConfig, Stage, StageRiskConfig,
StageStateConfig,
};
let n_blks = block_durations_hours.len();
assert!(n_blks > 0, "must have at least one block");
let bus = Bus {
id: EntityId(1),
name: "B1".to_string(),
deficit_segments: vec![DeficitSegment {
depth_mw: None,
cost_per_mwh: 500.0,
}],
excess_cost: 0.0,
};
let hydro = Hydro {
id: EntityId(2),
name: "FPHA1".to_string(),
bus_id: EntityId(1),
downstream_id: None,
entry_stage_id: None,
exit_stage_id: None,
min_storage_hm3: 0.0,
max_storage_hm3: 500.0,
min_outflow_m3s: 0.0,
max_outflow_m3s: None,
generation_model: HydroGenerationModel::Fpha,
min_turbined_m3s: 0.0,
max_turbined_m3s: 150.0,
min_generation_mw: 0.0,
max_generation_mw: 300.0,
tailrace: None,
hydraulic_losses: None,
efficiency: None,
evaporation_coefficients_mm: None,
evaporation_reference_volumes_hm3: None,
diversion: None,
filling: None,
penalties: HydroPenalties {
spillage_cost: 0.01,
diversion_cost: 0.0,
fpha_turbined_cost,
storage_violation_below_cost: 0.0,
filling_target_violation_cost: 0.0,
turbined_violation_below_cost: 0.0,
outflow_violation_below_cost: 0.0,
outflow_violation_above_cost: 0.0,
generation_violation_below_cost: 0.0,
evaporation_violation_cost: 0.0,
water_withdrawal_violation_cost: 0.0,
},
};
let blocks: Vec<Block> = block_durations_hours
.iter()
.enumerate()
.map(|(i, &hours)| Block {
index: i,
name: format!("BLK{i}"),
duration_hours: hours,
})
.collect();
let stages: Vec<Stage> = vec![Stage {
index: 0,
id: 0,
start_date: NaiveDate::from_ymd_opt(2024, 1, 1).unwrap(),
end_date: NaiveDate::from_ymd_opt(2024, 2, 1).unwrap(),
season_id: None,
blocks,
block_mode: BlockMode::Parallel,
state_config: StageStateConfig {
storage: true,
inflow_lags: false,
},
risk_config: StageRiskConfig::Expectation,
scenario_config: ScenarioSourceConfig {
branching_factor: 1,
noise_method: NoiseMethod::Saa,
},
}];
let inflow_models: Vec<InflowModel> = vec![InflowModel {
hydro_id: EntityId(2),
stage_id: 0,
mean_m3s: 80.0,
std_m3s: 20.0,
ar_coefficients: vec![],
residual_std_ratio: 1.0,
}];
let load_models: Vec<LoadModel> = vec![LoadModel {
bus_id: EntityId(1),
stage_id: 0,
mean_mw: 200.0,
std_mw: 0.0,
}];
let bounds = ResolvedBounds::new(
1,
0,
0,
0,
0,
1,
HydroStageBounds {
min_storage_hm3: 0.0,
max_storage_hm3: 500.0,
min_turbined_m3s: 0.0,
max_turbined_m3s: 150.0,
min_outflow_m3s: 0.0,
max_outflow_m3s: None,
min_generation_mw: 0.0,
max_generation_mw: 300.0,
max_diversion_m3s: None,
filling_inflow_m3s: 0.0,
water_withdrawal_m3s: 0.0,
},
ThermalStageBounds {
min_generation_mw: 0.0,
max_generation_mw: 0.0,
},
LineStageBounds {
direct_mw: 0.0,
reverse_mw: 0.0,
},
PumpingStageBounds {
min_flow_m3s: 0.0,
max_flow_m3s: 0.0,
},
ContractStageBounds {
min_mw: 0.0,
max_mw: 0.0,
price_per_mwh: 0.0,
},
);
let penalties = ResolvedPenalties::new(
1,
1,
0,
0,
1,
HydroStagePenalties {
spillage_cost: 0.01,
diversion_cost: 0.0,
fpha_turbined_cost,
storage_violation_below_cost: 0.0,
filling_target_violation_cost: 0.0,
turbined_violation_below_cost: 0.0,
outflow_violation_below_cost: 0.0,
outflow_violation_above_cost: 0.0,
generation_violation_below_cost: 0.0,
evaporation_violation_cost: 0.0,
water_withdrawal_violation_cost: 0.0,
},
BusStagePenalties { excess_cost: 0.0 },
LineStagePenalties { exchange_cost: 0.0 },
NcsStagePenalties {
curtailment_cost: 0.0,
},
);
let system = SystemBuilder::new()
.buses(vec![bus])
.hydros(vec![hydro])
.stages(stages)
.inflow_models(inflow_models)
.load_models(load_models)
.bounds(bounds)
.penalties(penalties)
.build()
.expect("fpha_system_with_turbined_cost: valid");
let plane = FphaPlane {
intercept: 10.0,
gamma_v: 0.5,
gamma_q: 2.0,
gamma_s: 0.1,
};
let planes = vec![plane; n_planes];
let models = vec![vec![ResolvedProductionModel::Fpha {
planes,
turbined_cost: 0.0,
}]];
let production = ProductionModelSet::new(models, 1, 1);
(system, production)
}
#[test]
fn fpha_turbined_cost_applied_to_fpha_turbine_column() {
let (system, production) = fpha_system_with_turbined_cost(3, 0.5, &[720.0]);
let result = build_stage_templates(
&system,
&no_penalty_config(),
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&production,
&default_evaporation(&system),
)
.expect("fpha system builds ok");
let t = &result.templates[0];
let turbine_col = 4_usize;
let expected = 0.5 * 720.0 / COST_SCALE_FACTOR;
assert!(
(t.objective[turbine_col] - expected).abs() < 1e-15,
"FPHA turbine col objective: expected {expected}, got {}",
t.objective[turbine_col]
);
}
#[test]
fn constant_hydro_turbine_column_has_zero_objective() {
let system = one_hydro_system(1, 0);
let result = build_stage_templates(
&system,
&no_penalty_config(),
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&default_production(&system),
&default_evaporation(&system),
)
.expect("constant productivity ok");
let t = &result.templates[0];
let turbine_col = 4_usize;
assert_eq!(
t.objective[turbine_col], 0.0,
"constant hydro turbine column must have zero objective, got {}",
t.objective[turbine_col]
);
}
#[test]
fn fpha_turbined_cost_multi_block_uses_per_block_hours() {
let (system, production) = fpha_system_with_turbined_cost(3, 1.0, &[300.0, 420.0]);
let result = build_stage_templates(
&system,
&no_penalty_config(),
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&production,
&default_evaporation(&system),
)
.expect("fpha multi-block system builds ok");
let t = &result.templates[0];
let col_blk0 = 4_usize;
let col_blk1 = 5_usize;
assert!(
(t.objective[col_blk0] - 300.0 / COST_SCALE_FACTOR).abs() < 1e-15,
"block 0 turbine objective: expected {}, got {}",
300.0 / COST_SCALE_FACTOR,
t.objective[col_blk0]
);
assert!(
(t.objective[col_blk1] - 420.0 / COST_SCALE_FACTOR).abs() < 1e-15,
"block 1 turbine objective: expected {}, got {}",
420.0 / COST_SCALE_FACTOR,
t.objective[col_blk1]
);
}
#[test]
fn fpha_turbined_cost_mixed_system_only_fpha_hydros_carry_cost() {
let (system, production) = four_hydro_mixed_system();
let hydro_pen = HydroStagePenalties {
spillage_cost: 0.01,
diversion_cost: 0.0,
fpha_turbined_cost: 1.0,
storage_violation_below_cost: 0.0,
filling_target_violation_cost: 0.0,
turbined_violation_below_cost: 0.0,
outflow_violation_below_cost: 0.0,
outflow_violation_above_cost: 0.0,
generation_violation_below_cost: 0.0,
evaporation_violation_cost: 0.0,
water_withdrawal_violation_cost: 0.0,
};
let penalties = ResolvedPenalties::new(
4,
1,
0,
0,
1,
hydro_pen,
BusStagePenalties { excess_cost: 0.0 },
LineStagePenalties { exchange_cost: 0.0 },
NcsStagePenalties {
curtailment_cost: 0.0,
},
);
let system = SystemBuilder::new()
.buses(system.buses().to_vec())
.hydros(system.hydros().to_vec())
.stages(system.stages().to_vec())
.inflow_models(system.inflow_models().to_vec())
.load_models(system.load_models().to_vec())
.bounds(system.bounds().clone())
.penalties(penalties)
.build()
.expect("mixed system with turbined cost");
let result = build_stage_templates(
&system,
&no_penalty_config(),
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&production,
&default_evaporation(&system),
)
.expect("mixed system builds");
let t = &result.templates[0];
let col_turbine_start = 13;
let block_hours = 744.0;
assert!(
t.objective[col_turbine_start].abs() < 1e-12,
"constant hydro 0 turbine objective should be 0.0, got {}",
t.objective[col_turbine_start]
);
assert!(
t.objective[col_turbine_start + 1].abs() < 1e-12,
"constant hydro 1 turbine objective should be 0.0, got {}",
t.objective[col_turbine_start + 1]
);
let expected_fpha = block_hours / COST_SCALE_FACTOR;
assert!(
(t.objective[col_turbine_start + 2] - expected_fpha).abs() < 1e-15,
"FPHA hydro 2 turbine objective should be {expected_fpha}, got {}",
t.objective[col_turbine_start + 2]
);
assert!(
(t.objective[col_turbine_start + 3] - expected_fpha).abs() < 1e-15,
"FPHA hydro 3 turbine objective should be {expected_fpha}, got {}",
t.objective[col_turbine_start + 3]
);
}
#[test]
fn load_balance_rhs_matches_load_model_mean_mw() {
let system = one_bus_system(1);
let result = build_stage_templates(
&system,
&no_penalty_config(),
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&default_production(&system),
&default_evaporation(&system),
)
.expect("constant productivity ok");
let t = &result.templates[0];
let load_row = 0;
assert_eq!(
t.row_lower[load_row], 100.0,
"row_lower for load balance must be mean_mw"
);
assert_eq!(
t.row_upper[load_row], 100.0,
"row_upper for load balance must be mean_mw"
);
}
#[test]
fn multiple_stages_produce_same_count_templates_and_base_rows() {
let system = one_hydro_system(3, 1);
let result = build_stage_templates(
&system,
&no_penalty_config(),
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&default_production(&system),
&default_evaporation(&system),
)
.expect("constant productivity ok");
assert_eq!(result.templates.len(), 3);
assert_eq!(result.base_rows.len(), 3);
}
#[test]
fn stage_templates_clone_and_debug() {
let system = one_hydro_system(1, 0);
let result = build_stage_templates(
&system,
&no_penalty_config(),
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&default_production(&system),
&default_evaporation(&system),
)
.expect("constant productivity ok");
let cloned = result.clone();
assert_eq!(cloned.templates.len(), result.templates.len());
let s = format!("{result:?}");
assert!(s.contains("StageTemplates"));
}
#[test]
#[allow(clippy::too_many_lines)]
fn test_fpha_model_accepted() {
use chrono::NaiveDate;
use cobre_core::entities::hydro::{Hydro, HydroGenerationModel, HydroPenalties};
use cobre_core::scenario::{InflowModel, LoadModel};
use cobre_core::temporal::{
Block, BlockMode, NoiseMethod, ScenarioSourceConfig, Stage, StageRiskConfig,
StageStateConfig,
};
let bus = Bus {
id: EntityId(1),
name: "B1".to_string(),
deficit_segments: vec![DeficitSegment {
depth_mw: None,
cost_per_mwh: 500.0,
}],
excess_cost: 0.0,
};
let hydro = Hydro {
id: EntityId(5),
name: "Tucurui".to_string(),
bus_id: EntityId(1),
downstream_id: None,
entry_stage_id: None,
exit_stage_id: None,
min_storage_hm3: 0.0,
max_storage_hm3: 200.0,
min_outflow_m3s: 0.0,
max_outflow_m3s: None,
generation_model: HydroGenerationModel::Fpha,
min_turbined_m3s: 0.0,
max_turbined_m3s: 100.0,
min_generation_mw: 0.0,
max_generation_mw: 250.0,
tailrace: None,
hydraulic_losses: None,
efficiency: None,
evaporation_coefficients_mm: None,
evaporation_reference_volumes_hm3: None,
diversion: None,
filling: None,
penalties: HydroPenalties {
spillage_cost: 0.01,
diversion_cost: 0.0,
fpha_turbined_cost: 0.0,
storage_violation_below_cost: 0.0,
filling_target_violation_cost: 0.0,
turbined_violation_below_cost: 0.0,
outflow_violation_below_cost: 0.0,
outflow_violation_above_cost: 0.0,
generation_violation_below_cost: 0.0,
evaporation_violation_cost: 0.0,
water_withdrawal_violation_cost: 0.0,
},
};
let stages: Vec<Stage> = vec![Stage {
index: 0,
id: 0,
start_date: NaiveDate::from_ymd_opt(2024, 1, 1).expect("valid date"),
end_date: NaiveDate::from_ymd_opt(2024, 2, 1).expect("valid date"),
season_id: None,
blocks: vec![Block {
index: 0,
name: "S".to_string(),
duration_hours: 744.0,
}],
block_mode: BlockMode::Parallel,
state_config: StageStateConfig {
storage: true,
inflow_lags: false,
},
risk_config: StageRiskConfig::Expectation,
scenario_config: ScenarioSourceConfig {
branching_factor: 1,
noise_method: NoiseMethod::Saa,
},
}];
let inflow_models: Vec<InflowModel> = vec![InflowModel {
hydro_id: EntityId(5),
stage_id: 0,
mean_m3s: 80.0,
std_m3s: 20.0,
ar_coefficients: vec![],
residual_std_ratio: 1.0,
}];
let load_models: Vec<LoadModel> = vec![LoadModel {
bus_id: EntityId(1),
stage_id: 0,
mean_mw: 100.0,
std_mw: 0.0,
}];
let bounds = ResolvedBounds::new(
1,
0,
0,
0,
0,
1,
default_hydro_bounds(),
ThermalStageBounds {
min_generation_mw: 0.0,
max_generation_mw: 0.0,
},
LineStageBounds {
direct_mw: 0.0,
reverse_mw: 0.0,
},
PumpingStageBounds {
min_flow_m3s: 0.0,
max_flow_m3s: 0.0,
},
ContractStageBounds {
min_mw: 0.0,
max_mw: 0.0,
price_per_mwh: 0.0,
},
);
let penalties = ResolvedPenalties::new(
1,
1,
0,
0,
1,
default_hydro_penalties(),
BusStagePenalties { excess_cost: 0.0 },
LineStagePenalties { exchange_cost: 0.0 },
NcsStagePenalties {
curtailment_cost: 0.0,
},
);
let system = cobre_core::SystemBuilder::new()
.buses(vec![bus])
.hydros(vec![hydro])
.stages(stages)
.inflow_models(inflow_models)
.load_models(load_models)
.bounds(bounds)
.penalties(penalties)
.build()
.expect("test_fpha_model_accepted: valid system");
let production = PrepareHydroModelsResult::default_from_system(&system).production;
let result = build_stage_templates(
&system,
&no_penalty_config(),
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&production,
&default_evaporation(&system),
);
assert!(
result.is_ok(),
"Fpha entity model with ConstantProductivity resolved model must now succeed: {result:?}"
);
if let Err(ref e) = result {
let msg = e.to_string();
assert!(
!msg.contains("Tucurui"),
"unexpected error for Tucurui: {msg}"
);
}
}
#[test]
fn test_constant_productivity_accepted() {
let system = one_hydro_system(1, 0);
let result = build_stage_templates(
&system,
&no_penalty_config(),
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&default_production(&system),
&default_evaporation(&system),
);
assert!(
result.is_ok(),
"ConstantProductivity system must return Ok, got: {result:?}"
);
assert_eq!(
result.expect("accepted").templates.len(),
1,
"one study stage → one template"
);
}
#[test]
fn test_penalty_columns_added() {
let system = one_hydro_system(1, 0);
let without = build_stage_templates(
&system,
&no_penalty_config(),
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&default_production(&system),
&default_evaporation(&system),
)
.expect("constant productivity ok");
let with_p = build_stage_templates(
&system,
&penalty_config(1000.0),
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&default_production(&system),
&default_evaporation(&system),
)
.expect("constant productivity ok");
assert_eq!(
with_p.templates[0].num_cols,
without.templates[0].num_cols + 1,
"penalty method must add exactly n_hydros extra columns"
);
}
#[test]
fn test_penalty_columns_added_3_hydros() {
let system = one_bus_system(1);
let with_p = build_stage_templates(
&system,
&penalty_config(1000.0),
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&default_production(&system),
&default_evaporation(&system),
)
.expect("constant productivity ok");
let without = build_stage_templates(
&system,
&no_penalty_config(),
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&default_production(&system),
&default_evaporation(&system),
)
.expect("constant productivity ok");
assert_eq!(
with_p.templates[0].num_cols, without.templates[0].num_cols,
"no slack columns when n_hydros == 0, even with penalty config"
);
}
#[test]
fn test_penalty_objective_coefficient() {
let system = one_hydro_system(1, 0);
let config = penalty_config(1000.0);
let result = build_stage_templates(
&system,
&config,
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&default_production(&system),
&default_evaporation(&system),
)
.expect("constant productivity ok");
let t = &result.templates[0];
let slack_col = t.num_cols - 1 - t.n_hydro; let expected_obj = 1000.0 * 744.0 / COST_SCALE_FACTOR;
assert!(
(t.objective[slack_col] - expected_obj).abs() < 1e-12,
"expected objective {expected_obj}, got {}",
t.objective[slack_col]
);
}
#[test]
fn test_no_penalty_columns_when_none() {
let system = one_hydro_system(1, 2);
let result = build_stage_templates(
&system,
&no_penalty_config(),
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&default_production(&system),
&default_evaporation(&system),
)
.expect("constant productivity ok");
let t = &result.templates[0];
assert_eq!(
t.num_cols, 12,
"method=none must not add extra penalty columns"
);
assert_eq!(t.num_rows, 6, "method=none must not add extra penalty rows");
}
#[test]
#[allow(clippy::cast_sign_loss)]
fn test_penalty_slack_in_water_balance() {
let system = one_hydro_system(1, 0);
let config = penalty_config(1000.0);
let result = build_stage_templates(
&system,
&config,
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&default_production(&system),
&default_evaporation(&system),
)
.expect("constant productivity ok");
let t = &result.templates[0];
let slack_col = t.num_cols - 1 - t.n_hydro;
let water_balance_row = 2_usize;
let col_start = t.col_starts[slack_col] as usize;
let col_end = t.col_starts[slack_col + 1] as usize;
let found = t.row_indices[col_start..col_end]
.iter()
.zip(&t.values[col_start..col_end])
.any(|(&r, &v)| r as usize == water_balance_row && v.abs() > 1e-12);
assert!(
found,
"slack column must have a non-zero entry in the water balance row"
);
}
#[test]
fn test_penalty_slack_bounds() {
let system = one_hydro_system(1, 0);
let config = penalty_config(1000.0);
let result = build_stage_templates(
&system,
&config,
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&default_production(&system),
&default_evaporation(&system),
)
.expect("constant productivity ok");
let t = &result.templates[0];
let slack_col = t.num_cols - 1 - t.n_hydro;
assert_eq!(t.col_lower[slack_col], 0.0, "slack lower bound must be 0.0");
assert!(
t.col_upper[slack_col].is_infinite() && t.col_upper[slack_col] > 0.0,
"slack upper bound must be +infinity"
);
}
#[test]
#[allow(clippy::cast_sign_loss)]
fn test_penalty_water_balance_coefficient_value() {
let system = one_hydro_system(1, 0);
let config = penalty_config(1000.0);
let result = build_stage_templates(
&system,
&config,
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&default_production(&system),
&default_evaporation(&system),
)
.expect("constant productivity ok");
let t = &result.templates[0];
let slack_col = t.num_cols - 1 - t.n_hydro;
let water_balance_row = 2_usize; let zeta = 744.0 * (3_600.0 / 1_000_000.0);
let expected_coeff = -zeta;
let col_start = t.col_starts[slack_col] as usize;
let col_end = t.col_starts[slack_col + 1] as usize;
let coeff = t.row_indices[col_start..col_end]
.iter()
.zip(&t.values[col_start..col_end])
.find(|&(&r, _)| r as usize == water_balance_row)
.map(|(_, &v)| v);
assert!(
coeff.is_some(),
"slack column must have an entry in the water balance row"
);
let coeff = coeff.unwrap();
assert!(
(coeff - expected_coeff).abs() < 1e-9,
"expected coefficient {expected_coeff:.9}, got {coeff:.9}"
);
}
#[test]
fn test_penalty_multi_stage_consistent() {
let system = one_hydro_system(3, 1);
let config = penalty_config(2000.0);
let result = build_stage_templates(
&system,
&config,
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&default_production(&system),
&default_evaporation(&system),
)
.expect("constant productivity ok");
assert_eq!(result.templates.len(), 3);
let base_cols = result.templates[0].num_cols;
for t in &result.templates {
assert_eq!(
t.num_cols, base_cols,
"all stages must have the same column count"
);
}
}
#[test]
fn test_penalty_slack_absorbs_negative_inflow() {
use cobre_solver::{HighsSolver, RowBatch, SolverInterface};
let system = one_hydro_system(1, 0);
let config = penalty_config(1000.0);
let result = build_stage_templates(
&system,
&config,
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&default_production(&system),
&default_evaporation(&system),
)
.expect("constant productivity ok");
let template = &result.templates[0];
let col_inflow_slack_start = template.num_cols - 1 - template.n_hydro;
let base_row = result.base_rows[0];
let water_balance_row = base_row;
let mut solver = HighsSolver::new().expect("HighsSolver::new must succeed");
solver.load_model(template);
let empty_cuts = RowBatch {
num_rows: 0,
row_starts: vec![0_i32],
col_indices: vec![],
values: vec![],
row_lower: vec![],
row_upper: vec![],
};
solver.add_rows(&empty_cuts);
let initial_storage = 100.0_f64;
let negative_noise = -5.0_f64;
solver.set_row_bounds(
&[0, water_balance_row],
&[initial_storage, negative_noise],
&[initial_storage, negative_noise],
);
let view = solver
.solve()
.expect("LP must be feasible with inflow slack active");
let primal = view.primal;
assert!(
primal[col_inflow_slack_start] > 0.0,
"inflow slack must be positive when noise is negative, got {}",
primal[col_inflow_slack_start]
);
assert!(
view.objective > 0.0,
"objective must include a positive penalty contribution, got {}",
view.objective
);
}
#[allow(
clippy::cast_possible_truncation,
clippy::cast_possible_wrap,
clippy::too_many_lines
)]
fn two_bus_system_with_stochastic_load(
n_stages: usize,
n_hydros_in_system: usize,
n_blocks: usize,
) -> cobre_core::System {
use chrono::NaiveDate;
use cobre_core::entities::hydro::{Hydro, HydroGenerationModel, HydroPenalties};
use cobre_core::scenario::{InflowModel, LoadModel};
use cobre_core::temporal::{
Block, BlockMode, NoiseMethod, ScenarioSourceConfig, Stage, StageRiskConfig,
StageStateConfig,
};
let bus1 = Bus {
id: EntityId(10),
name: "B1".to_string(),
deficit_segments: vec![DeficitSegment {
depth_mw: None,
cost_per_mwh: 500.0,
}],
excess_cost: 0.0,
};
let bus2 = Bus {
id: EntityId(20),
name: "B2".to_string(),
deficit_segments: vec![DeficitSegment {
depth_mw: None,
cost_per_mwh: 500.0,
}],
excess_cost: 0.0,
};
let blocks: Vec<_> = (0..n_blocks)
.map(|b| Block {
index: b,
name: format!("B{b}"),
duration_hours: 240.0,
})
.collect();
let stages: Vec<Stage> = (0..n_stages)
.map(|i| Stage {
index: i,
id: i as i32,
start_date: NaiveDate::from_ymd_opt(2024, 1, 1).unwrap(),
end_date: NaiveDate::from_ymd_opt(2024, 2, 1).unwrap(),
season_id: None,
blocks: blocks.clone(),
block_mode: BlockMode::Parallel,
state_config: StageStateConfig {
storage: false,
inflow_lags: false,
},
risk_config: StageRiskConfig::Expectation,
scenario_config: ScenarioSourceConfig {
branching_factor: 1,
noise_method: NoiseMethod::Saa,
},
})
.collect();
let hydros: Vec<Hydro> = (0..n_hydros_in_system)
.map(|h| Hydro {
id: EntityId((h + 100) as i32),
name: format!("H{h}"),
bus_id: EntityId(10),
downstream_id: None,
entry_stage_id: None,
exit_stage_id: None,
min_storage_hm3: 0.0,
max_storage_hm3: 200.0,
min_outflow_m3s: 0.0,
max_outflow_m3s: None,
generation_model: HydroGenerationModel::ConstantProductivity {
productivity_mw_per_m3s: 1.0,
},
min_turbined_m3s: 0.0,
max_turbined_m3s: 50.0,
min_generation_mw: 0.0,
max_generation_mw: 50.0,
tailrace: None,
hydraulic_losses: None,
efficiency: None,
evaporation_coefficients_mm: None,
evaporation_reference_volumes_hm3: None,
diversion: None,
filling: None,
penalties: HydroPenalties {
spillage_cost: 0.0,
diversion_cost: 0.0,
fpha_turbined_cost: 0.0,
storage_violation_below_cost: 0.0,
filling_target_violation_cost: 0.0,
turbined_violation_below_cost: 0.0,
outflow_violation_below_cost: 0.0,
outflow_violation_above_cost: 0.0,
generation_violation_below_cost: 0.0,
evaporation_violation_cost: 0.0,
water_withdrawal_violation_cost: 0.0,
},
})
.collect();
let inflow_models: Vec<InflowModel> = hydros
.iter()
.flat_map(|h| {
(0..n_stages).map(move |s| InflowModel {
hydro_id: h.id,
stage_id: s as i32,
mean_m3s: 50.0,
std_m3s: 10.0,
ar_coefficients: vec![],
residual_std_ratio: 1.0,
})
})
.collect();
let load_models: Vec<LoadModel> = (0..n_stages)
.flat_map(|s| {
[
LoadModel {
bus_id: EntityId(10),
stage_id: s as i32,
mean_mw: 80.0,
std_mw: 0.0, },
LoadModel {
bus_id: EntityId(20),
stage_id: s as i32,
mean_mw: 120.0,
std_mw: 15.0, },
]
})
.collect();
let n_st = n_stages.max(1);
let bounds = ResolvedBounds::new(
n_hydros_in_system,
0,
0,
0,
0,
n_st,
default_hydro_bounds(),
ThermalStageBounds {
min_generation_mw: 0.0,
max_generation_mw: 0.0,
},
LineStageBounds {
direct_mw: 0.0,
reverse_mw: 0.0,
},
PumpingStageBounds {
min_flow_m3s: 0.0,
max_flow_m3s: 0.0,
},
ContractStageBounds {
min_mw: 0.0,
max_mw: 0.0,
price_per_mwh: 0.0,
},
);
let penalties = ResolvedPenalties::new(
n_hydros_in_system,
2,
0,
0,
n_st,
default_hydro_penalties(),
BusStagePenalties { excess_cost: 0.0 },
LineStagePenalties { exchange_cost: 0.0 },
NcsStagePenalties {
curtailment_cost: 0.0,
},
);
let mut builder = SystemBuilder::new()
.buses(vec![bus1, bus2])
.stages(stages)
.load_models(load_models)
.bounds(bounds)
.penalties(penalties);
if !hydros.is_empty() {
builder = builder.hydros(hydros).inflow_models(inflow_models);
}
builder.build().expect("two_bus_system: valid")
}
#[test]
fn stage_templates_load_balance_row_starts_correct() {
let system = two_bus_system_with_stochastic_load(2, 2, 3);
let result = build_stage_templates(
&system,
&no_penalty_config(),
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&default_production(&system),
&default_evaporation(&system),
)
.expect("constant productivity ok");
assert_eq!(
result.load_balance_row_starts.len(),
result.templates.len(),
"load_balance_row_starts length must match templates length"
);
let expected_row_start = result.base_rows[0] + 2; assert_eq!(
result.load_balance_row_starts[0], expected_row_start,
"load_balance_row_starts[0] must equal row_water_balance_start + n_hydros"
);
assert_eq!(
result.load_balance_row_starts[0], result.load_balance_row_starts[1],
"identical stages share the same load balance row start"
);
}
#[test]
fn stage_templates_n_load_buses_matches_stochastic_buses() {
let system = two_bus_system_with_stochastic_load(1, 0, 1);
let result = build_stage_templates(
&system,
&no_penalty_config(),
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&default_production(&system),
&default_evaporation(&system),
)
.expect("constant productivity ok");
assert_eq!(
result.n_load_buses, 1,
"only B2 has std_mw > 0 → n_load_buses must be 1"
);
assert_eq!(
result.load_bus_indices.len(),
1,
"load_bus_indices must have exactly one entry"
);
assert_eq!(
result.load_bus_indices[0], 1,
"B2 is at buses slice index 1 (buses are [B1(10), B2(20)])"
);
}
#[test]
fn stage_templates_no_load_buses_gives_zero() {
let system = one_bus_system(2);
let result = build_stage_templates(
&system,
&no_penalty_config(),
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&default_production(&system),
&default_evaporation(&system),
)
.expect("constant productivity ok");
assert_eq!(
result.n_load_buses, 0,
"system with std_mw = 0 everywhere must give n_load_buses = 0"
);
assert!(
result.load_bus_indices.is_empty(),
"load_bus_indices must be empty when n_load_buses = 0"
);
assert_eq!(
result.load_balance_row_starts.len(),
result.templates.len(),
"load_balance_row_starts length must always match templates length"
);
}
#[allow(clippy::cast_sign_loss)] fn csc_entry(tmpl: &cobre_solver::StageTemplate, col: usize, row: usize) -> Option<f64> {
let start = tmpl.col_starts[col] as usize;
let end = tmpl.col_starts[col + 1] as usize;
for pos in start..end {
if tmpl.row_indices[pos] as usize == row {
return Some(tmpl.values[pos]);
}
}
None
}
#[allow(
clippy::cast_possible_truncation,
clippy::cast_possible_wrap,
clippy::too_many_lines
)]
fn one_fpha_hydro_system(n_planes: usize) -> (cobre_core::System, ProductionModelSet) {
use chrono::NaiveDate;
use cobre_core::entities::hydro::{Hydro, HydroGenerationModel, HydroPenalties};
use cobre_core::scenario::{InflowModel, LoadModel};
use cobre_core::temporal::{
Block, BlockMode, NoiseMethod, ScenarioSourceConfig, Stage, StageRiskConfig,
StageStateConfig,
};
let bus = Bus {
id: EntityId(1),
name: "B1".to_string(),
deficit_segments: vec![DeficitSegment {
depth_mw: None,
cost_per_mwh: 500.0,
}],
excess_cost: 0.0,
};
let hydro = Hydro {
id: EntityId(2),
name: "FPHA1".to_string(),
bus_id: EntityId(1),
downstream_id: None,
entry_stage_id: None,
exit_stage_id: None,
min_storage_hm3: 0.0,
max_storage_hm3: 500.0,
min_outflow_m3s: 0.0,
max_outflow_m3s: None,
generation_model: HydroGenerationModel::Fpha,
min_turbined_m3s: 0.0,
max_turbined_m3s: 150.0,
min_generation_mw: 0.0,
max_generation_mw: 300.0,
tailrace: None,
hydraulic_losses: None,
efficiency: None,
evaporation_coefficients_mm: None,
evaporation_reference_volumes_hm3: None,
diversion: None,
filling: None,
penalties: HydroPenalties {
spillage_cost: 0.01,
diversion_cost: 0.0,
fpha_turbined_cost: 0.0,
storage_violation_below_cost: 0.0,
filling_target_violation_cost: 0.0,
turbined_violation_below_cost: 0.0,
outflow_violation_below_cost: 0.0,
outflow_violation_above_cost: 0.0,
generation_violation_below_cost: 0.0,
evaporation_violation_cost: 0.0,
water_withdrawal_violation_cost: 0.0,
},
};
let stages: Vec<Stage> = vec![Stage {
index: 0,
id: 0,
start_date: NaiveDate::from_ymd_opt(2024, 1, 1).unwrap(),
end_date: NaiveDate::from_ymd_opt(2024, 2, 1).unwrap(),
season_id: None,
blocks: vec![Block {
index: 0,
name: "S".to_string(),
duration_hours: 744.0,
}],
block_mode: BlockMode::Parallel,
state_config: StageStateConfig {
storage: true,
inflow_lags: false,
},
risk_config: StageRiskConfig::Expectation,
scenario_config: ScenarioSourceConfig {
branching_factor: 1,
noise_method: NoiseMethod::Saa,
},
}];
let inflow_models: Vec<InflowModel> = vec![InflowModel {
hydro_id: EntityId(2),
stage_id: 0,
mean_m3s: 80.0,
std_m3s: 20.0,
ar_coefficients: vec![],
residual_std_ratio: 1.0,
}];
let load_models: Vec<LoadModel> = vec![LoadModel {
bus_id: EntityId(1),
stage_id: 0,
mean_mw: 200.0,
std_mw: 0.0,
}];
let bounds = ResolvedBounds::new(
1,
0,
0,
0,
0,
1,
HydroStageBounds {
min_storage_hm3: 0.0,
max_storage_hm3: 500.0,
min_turbined_m3s: 0.0,
max_turbined_m3s: 150.0,
min_outflow_m3s: 0.0,
max_outflow_m3s: None,
min_generation_mw: 0.0,
max_generation_mw: 300.0,
max_diversion_m3s: None,
filling_inflow_m3s: 0.0,
water_withdrawal_m3s: 0.0,
},
ThermalStageBounds {
min_generation_mw: 0.0,
max_generation_mw: 0.0,
},
LineStageBounds {
direct_mw: 0.0,
reverse_mw: 0.0,
},
PumpingStageBounds {
min_flow_m3s: 0.0,
max_flow_m3s: 0.0,
},
ContractStageBounds {
min_mw: 0.0,
max_mw: 0.0,
price_per_mwh: 0.0,
},
);
let penalties = ResolvedPenalties::new(
1,
1,
0,
0,
1,
default_hydro_penalties(),
BusStagePenalties { excess_cost: 0.0 },
LineStagePenalties { exchange_cost: 0.0 },
NcsStagePenalties {
curtailment_cost: 0.0,
},
);
let system = SystemBuilder::new()
.buses(vec![bus])
.hydros(vec![hydro])
.stages(stages)
.inflow_models(inflow_models)
.load_models(load_models)
.bounds(bounds)
.penalties(penalties)
.build()
.expect("one_fpha_hydro_system: valid");
let plane = FphaPlane {
intercept: 10.0,
gamma_v: 0.5,
gamma_q: 2.0,
gamma_s: 0.1,
};
let planes = vec![plane; n_planes];
let models = vec![vec![ResolvedProductionModel::Fpha {
planes,
turbined_cost: 0.0,
}]];
let production = ProductionModelSet::new(models, 1, 1);
(system, production)
}
#[allow(
clippy::cast_possible_truncation,
clippy::cast_possible_wrap,
clippy::too_many_lines
)]
fn four_hydro_mixed_system() -> (cobre_core::System, ProductionModelSet) {
use chrono::NaiveDate;
use cobre_core::entities::hydro::{Hydro, HydroGenerationModel, HydroPenalties};
use cobre_core::scenario::{InflowModel, LoadModel};
use cobre_core::temporal::{
Block, BlockMode, NoiseMethod, ScenarioSourceConfig, Stage, StageRiskConfig,
StageStateConfig,
};
let bus = Bus {
id: EntityId(1),
name: "B1".to_string(),
deficit_segments: vec![DeficitSegment {
depth_mw: None,
cost_per_mwh: 500.0,
}],
excess_cost: 0.0,
};
let hydro_penalties = HydroPenalties {
spillage_cost: 0.01,
diversion_cost: 0.0,
fpha_turbined_cost: 0.0,
storage_violation_below_cost: 0.0,
filling_target_violation_cost: 0.0,
turbined_violation_below_cost: 0.0,
outflow_violation_below_cost: 0.0,
outflow_violation_above_cost: 0.0,
generation_violation_below_cost: 0.0,
evaporation_violation_cost: 0.0,
water_withdrawal_violation_cost: 0.0,
};
let make_hydro = |id: i32, gen_model: HydroGenerationModel| Hydro {
id: EntityId(id),
name: format!("H{id}"),
bus_id: EntityId(1),
downstream_id: None,
entry_stage_id: None,
exit_stage_id: None,
min_storage_hm3: 0.0,
max_storage_hm3: 200.0,
min_outflow_m3s: 0.0,
max_outflow_m3s: None,
generation_model: gen_model,
min_turbined_m3s: 0.0,
max_turbined_m3s: 100.0,
min_generation_mw: 0.0,
max_generation_mw: 250.0,
tailrace: None,
hydraulic_losses: None,
efficiency: None,
evaporation_coefficients_mm: None,
evaporation_reference_volumes_hm3: None,
diversion: None,
filling: None,
penalties: hydro_penalties,
};
let hydros = vec![
make_hydro(
100,
HydroGenerationModel::ConstantProductivity {
productivity_mw_per_m3s: 2.5,
},
),
make_hydro(
101,
HydroGenerationModel::ConstantProductivity {
productivity_mw_per_m3s: 3.0,
},
),
make_hydro(102, HydroGenerationModel::Fpha),
make_hydro(103, HydroGenerationModel::Fpha),
];
let stages: Vec<Stage> = vec![Stage {
index: 0,
id: 0,
start_date: NaiveDate::from_ymd_opt(2024, 1, 1).unwrap(),
end_date: NaiveDate::from_ymd_opt(2024, 2, 1).unwrap(),
season_id: None,
blocks: vec![Block {
index: 0,
name: "S".to_string(),
duration_hours: 744.0,
}],
block_mode: BlockMode::Parallel,
state_config: StageStateConfig {
storage: true,
inflow_lags: false,
},
risk_config: StageRiskConfig::Expectation,
scenario_config: ScenarioSourceConfig {
branching_factor: 1,
noise_method: NoiseMethod::Saa,
},
}];
let inflow_models: Vec<InflowModel> = hydros
.iter()
.map(|h| InflowModel {
hydro_id: h.id,
stage_id: 0,
mean_m3s: 80.0,
std_m3s: 20.0,
ar_coefficients: vec![],
residual_std_ratio: 1.0,
})
.collect();
let load_models: Vec<LoadModel> = vec![LoadModel {
bus_id: EntityId(1),
stage_id: 0,
mean_mw: 400.0,
std_mw: 0.0,
}];
let bounds = ResolvedBounds::new(
4,
0,
0,
0,
0,
1,
default_hydro_bounds(),
ThermalStageBounds {
min_generation_mw: 0.0,
max_generation_mw: 0.0,
},
LineStageBounds {
direct_mw: 0.0,
reverse_mw: 0.0,
},
PumpingStageBounds {
min_flow_m3s: 0.0,
max_flow_m3s: 0.0,
},
ContractStageBounds {
min_mw: 0.0,
max_mw: 0.0,
price_per_mwh: 0.0,
},
);
let penalties = ResolvedPenalties::new(
4,
1,
0,
0,
1,
default_hydro_penalties(),
BusStagePenalties { excess_cost: 0.0 },
LineStagePenalties { exchange_cost: 0.0 },
NcsStagePenalties {
curtailment_cost: 0.0,
},
);
let system = SystemBuilder::new()
.buses(vec![bus])
.hydros(hydros)
.stages(stages)
.inflow_models(inflow_models)
.load_models(load_models)
.bounds(bounds)
.penalties(penalties)
.build()
.expect("four_hydro_mixed_system: valid");
let plane = FphaPlane {
intercept: 10.0,
gamma_v: 0.5,
gamma_q: 2.0,
gamma_s: 0.1,
};
let fpha_planes = vec![plane; 3];
let models = vec![
vec![ResolvedProductionModel::ConstantProductivity { productivity: 2.5 }],
vec![ResolvedProductionModel::ConstantProductivity { productivity: 3.0 }],
vec![ResolvedProductionModel::Fpha {
planes: fpha_planes.clone(),
turbined_cost: 0.0,
}],
vec![ResolvedProductionModel::Fpha {
planes: fpha_planes,
turbined_cost: 0.0,
}],
];
let production = ProductionModelSet::new(models, 4, 1);
(system, production)
}
#[test]
fn fpha_ac1_dimensions_one_fpha_hydro_five_planes() {
let (system, production) = one_fpha_hydro_system(5);
let fpha_result = build_stage_templates(
&system,
&no_penalty_config(),
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&production,
&default_evaporation(&system),
)
.expect("FPHA system ok");
let const_result = build_stage_templates(
&system,
&no_penalty_config(),
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&default_production(&system),
&default_evaporation(&system),
)
.expect("constant productivity ok");
let fpha_tmpl = &fpha_result.templates[0];
let const_tmpl = &const_result.templates[0];
assert_eq!(
fpha_tmpl.num_cols,
const_tmpl.num_cols + 1,
"FPHA adds exactly 1 generation column (1 hydro * 1 block)"
);
assert_eq!(
fpha_tmpl.num_rows,
const_tmpl.num_rows + 5,
"FPHA adds exactly 5 constraint rows (5 planes * 1 block)"
);
}
#[test]
fn fpha_ac2_generation_column_entries() {
let n_planes = 5;
let (system, production) = one_fpha_hydro_system(n_planes);
let result = build_stage_templates(
&system,
&no_penalty_config(),
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&production,
&default_evaporation(&system),
)
.expect("FPHA system ok");
let tmpl = &result.templates[0];
let col_g = 9_usize;
let row_fpha_start = 4_usize;
for p in 0..n_planes {
let row = row_fpha_start + p;
let coeff = csc_entry(tmpl, col_g, row)
.unwrap_or_else(|| panic!("generation column missing entry in FPHA row {row}"));
assert!(
(coeff - 1.0).abs() < 1e-12,
"generation col FPHA row {row}: expected +1.0, got {coeff}"
);
}
let row_lb = 3_usize;
let lb_coeff = csc_entry(tmpl, col_g, row_lb).unwrap_or_else(|| {
panic!("generation column missing entry in load balance row {row_lb}")
});
assert!(
(lb_coeff - 1.0).abs() < 1e-12,
"generation col load balance row: expected +1.0, got {lb_coeff}"
);
}
#[test]
fn fpha_ac3_v_in_column_entries() {
let n_planes = 5;
let (system, production) = one_fpha_hydro_system(n_planes);
let result = build_stage_templates(
&system,
&no_penalty_config(),
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&production,
&default_evaporation(&system),
)
.expect("FPHA system ok");
let tmpl = &result.templates[0];
let col_v_in = 2_usize;
let row_fpha_start = 4_usize;
let expected = -0.5_f64 / 2.0;
for p in 0..n_planes {
let row = row_fpha_start + p;
let coeff = csc_entry(tmpl, col_v_in, row)
.unwrap_or_else(|| panic!("v_in column missing entry in FPHA row {row}"));
assert!(
(coeff - expected).abs() < 1e-12,
"v_in col FPHA row {row}: expected {expected}, got {coeff}"
);
}
}
#[test]
fn fpha_ac4_v_out_column_entries() {
let n_planes = 5;
let (system, production) = one_fpha_hydro_system(n_planes);
let result = build_stage_templates(
&system,
&no_penalty_config(),
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&production,
&default_evaporation(&system),
)
.expect("FPHA system ok");
let tmpl = &result.templates[0];
let col_v = 0_usize;
let row_fpha_start = 4_usize;
let expected = -0.5_f64 / 2.0;
for p in 0..n_planes {
let row = row_fpha_start + p;
let coeff = csc_entry(tmpl, col_v, row)
.unwrap_or_else(|| panic!("v column missing entry in FPHA row {row}"));
assert!(
(coeff - expected).abs() < 1e-12,
"v col FPHA row {row}: expected {expected}, got {coeff}"
);
}
}
#[test]
fn fpha_ac5_mixed_system_load_balance_uses_generation_col() {
let (system, production) = four_hydro_mixed_system();
let result = build_stage_templates(
&system,
&no_penalty_config(),
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&production,
&default_evaporation(&system),
)
.expect("mixed FPHA/constant system ok");
let tmpl = &result.templates[0];
assert_eq!(
tmpl.num_cols, 33,
"4-hydro mixed system: num_cols should be 33 (includes diversion and withdrawal slack)"
);
assert_eq!(
tmpl.num_rows, 19,
"4-hydro mixed system: num_rows should be 19 (includes z_inflow rows)"
);
let row_lb = 12_usize;
let col_g_fpha0 = 27_usize;
let g0_lb_coeff = csc_entry(tmpl, col_g_fpha0, row_lb).unwrap_or_else(|| {
panic!("FPHA hydro 0 generation column missing entry in load balance row {row_lb}")
});
assert!(
(g0_lb_coeff - 1.0).abs() < 1e-12,
"FPHA hydro 0 load balance: expected +1.0, got {g0_lb_coeff}"
);
let col_g_fpha1 = 28_usize;
let g1_lb_coeff = csc_entry(tmpl, col_g_fpha1, row_lb).unwrap_or_else(|| {
panic!("FPHA hydro 1 generation column missing entry in load balance row {row_lb}")
});
assert!(
(g1_lb_coeff - 1.0).abs() < 1e-12,
"FPHA hydro 1 load balance: expected +1.0, got {g1_lb_coeff}"
);
let col_turb_const = 13_usize;
let turb_lb_coeff = csc_entry(tmpl, col_turb_const, row_lb);
assert!(
turb_lb_coeff.is_some(),
"constant hydro 0 turbine col must appear in load balance row"
);
let expected_rho_coeff = 2.5_f64;
assert!(
(turb_lb_coeff.unwrap() - expected_rho_coeff).abs() < 1e-12,
"constant hydro 0 turbine: expected rho = {expected_rho_coeff}, got {turb_lb_coeff:?}"
);
}
fn fpha_solve_system() -> (cobre_core::System, ProductionModelSet) {
let planes = vec![
FphaPlane {
intercept: 300.0,
gamma_v: 1.0,
gamma_q: 3.0,
gamma_s: 0.0,
};
3
];
let models = vec![vec![ResolvedProductionModel::Fpha {
planes,
turbined_cost: 0.0,
}]];
let production = ProductionModelSet::new(models, 1, 1);
let (system, _) = one_fpha_hydro_system(3);
(system, production)
}
#[test]
fn fpha_solve_one_hydro_optimal() {
use cobre_solver::{HighsSolver, RowBatch, SolverInterface};
let (system, production) = fpha_solve_system();
let result = build_stage_templates(
&system,
&no_penalty_config(),
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&production,
&default_evaporation(&system),
)
.expect("FPHA template build must succeed");
let template = &result.templates[0];
let mut solver = HighsSolver::new().expect("HighsSolver::new must succeed");
solver.load_model(template);
let empty_cuts = RowBatch {
num_rows: 0,
row_starts: vec![0_i32],
col_indices: vec![],
values: vec![],
row_lower: vec![],
row_upper: vec![],
};
solver.add_rows(&empty_cuts);
let v_in = 100.0_f64;
solver.set_row_bounds(&[0], &[v_in], &[v_in]);
let view = solver
.solve()
.expect("FPHA LP must be feasible and optimal");
let col_g = 9_usize;
let generation = view.primal[col_g];
assert!(
generation > 0.0,
"FPHA generation must be strictly positive, got {generation}"
);
}
#[test]
fn fpha_solve_hyperplane_constraints_hold() {
use cobre_solver::{HighsSolver, RowBatch, SolverInterface};
let (system, production) = fpha_solve_system();
let planes = match production.model(0, 0) {
ResolvedProductionModel::Fpha { planes, .. } => planes.clone(),
ResolvedProductionModel::ConstantProductivity { .. } => {
panic!("expected Fpha model")
}
};
let result = build_stage_templates(
&system,
&no_penalty_config(),
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&production,
&default_evaporation(&system),
)
.expect("FPHA template build must succeed");
let template = &result.templates[0];
let mut solver = HighsSolver::new().expect("HighsSolver::new must succeed");
solver.load_model(template);
let empty_cuts = RowBatch {
num_rows: 0,
row_starts: vec![0_i32],
col_indices: vec![],
values: vec![],
row_lower: vec![],
row_upper: vec![],
};
solver.add_rows(&empty_cuts);
let v_in = 100.0_f64;
solver.set_row_bounds(&[0], &[v_in], &[v_in]);
let view = solver.solve().expect("FPHA LP must solve to optimal");
let primal = view.primal;
let col_v = 0_usize;
let col_v_in = 1_usize;
let col_q = 3_usize;
let col_s = 4_usize;
let col_g = 9_usize;
let g = primal[col_g];
let v = primal[col_v];
let v_in_sol = primal[col_v_in];
let q = primal[col_q];
let s = primal[col_s];
let v_avg = f64::midpoint(v, v_in_sol);
for (p_idx, plane) in planes.iter().enumerate() {
let rhs =
plane.intercept + plane.gamma_v * v_avg + plane.gamma_q * q + plane.gamma_s * s;
assert!(
g <= rhs + 1e-6,
"FPHA plane {p_idx}: g={g} must be <= rhs={rhs} \
(intercept={intercept}, gamma_v={gamma_v}, v_avg={v_avg}, \
gamma_q={gamma_q}, q={q}, gamma_s={gamma_s}, s={s})",
intercept = plane.intercept,
gamma_v = plane.gamma_v,
gamma_q = plane.gamma_q,
gamma_s = plane.gamma_s,
);
}
}
#[test]
fn fpha_solve_storage_fixing_dual_differs_from_constant() {
use cobre_solver::{HighsSolver, RowBatch, SolverInterface};
let (system, _) = one_fpha_hydro_system(1);
let tight_planes = vec![FphaPlane {
intercept: 0.0,
gamma_v: 0.5,
gamma_q: 1.0,
gamma_s: 0.0,
}];
let fpha_production = ProductionModelSet::new(
vec![vec![ResolvedProductionModel::Fpha {
planes: tight_planes,
turbined_cost: 0.0,
}]],
1,
1,
);
let fpha_result = build_stage_templates(
&system,
&no_penalty_config(),
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&fpha_production,
&default_evaporation(&system),
)
.expect("FPHA template build must succeed");
let const_production = default_production(&system);
let const_result = build_stage_templates(
&system,
&no_penalty_config(),
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&const_production,
&default_evaporation(&system),
)
.expect("constant productivity template build must succeed");
let solve_and_get_storage_dual = |template: &cobre_solver::StageTemplate| -> f64 {
let mut solver = HighsSolver::new().expect("HighsSolver::new must succeed");
solver.load_model(template);
let empty_cuts = RowBatch {
num_rows: 0,
row_starts: vec![0_i32],
col_indices: vec![],
values: vec![],
row_lower: vec![],
row_upper: vec![],
};
solver.add_rows(&empty_cuts);
let v_in = 100.0_f64;
solver.set_row_bounds(&[0], &[v_in], &[v_in]);
let view = solver.solve().expect("LP must solve to optimal");
view.dual[0]
};
let fpha_dual = solve_and_get_storage_dual(&fpha_result.templates[0]);
let const_dual = solve_and_get_storage_dual(&const_result.templates[0]);
assert!(
const_dual.abs() < 1e-6,
"constant-productivity dual must be ~0, got {const_dual}"
);
assert!(
fpha_dual.abs() > 1e-6,
"FPHA storage-fixing dual must be non-zero (FPHA v_in contribution \
must be present), got {fpha_dual}"
);
assert_ne!(
(fpha_dual * 1e6).round(),
(const_dual * 1e6).round(),
"storage-fixing dual must differ between FPHA ({fpha_dual}) and \
constant-productivity ({const_dual})"
);
}
#[test]
fn fpha_solve_mixed_system_optimal() {
use cobre_solver::{HighsSolver, RowBatch, SolverInterface};
let (system, production) = four_hydro_mixed_system();
let result = build_stage_templates(
&system,
&no_penalty_config(),
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&production,
&default_evaporation(&system),
)
.expect("mixed FPHA/constant system template build must succeed");
let template = &result.templates[0];
let mut solver = HighsSolver::new().expect("HighsSolver::new must succeed");
solver.load_model(template);
let empty_cuts = RowBatch {
num_rows: 0,
row_starts: vec![0_i32],
col_indices: vec![],
values: vec![],
row_lower: vec![],
row_upper: vec![],
};
solver.add_rows(&empty_cuts);
solver.set_row_bounds(
&[0, 1, 2, 3],
&[100.0, 100.0, 100.0, 100.0],
&[100.0, 100.0, 100.0, 100.0],
);
let view = solver
.solve()
.expect("mixed FPHA LP must be feasible and optimal");
assert!(
view.objective.is_finite(),
"objective must be finite, got {}",
view.objective
);
let col_g0 = 19_usize;
let col_g1 = 20_usize;
assert!(
view.primal[col_g0] >= 0.0,
"FPHA hydro 0 generation must be non-negative, got {}",
view.primal[col_g0]
);
assert!(
view.primal[col_g1] >= 0.0,
"FPHA hydro 1 generation must be non-negative, got {}",
view.primal[col_g1]
);
}
use cobre_solver::StageTemplate;
use crate::hydro_models::LinearizedEvaporation;
fn evap_set_for_system(
system: &cobre_core::System,
evap_indices: &[usize],
k_evap0_per_stage: &[f64],
) -> EvaporationModelSet {
let n_hydros = system.hydros().len();
let n_stages = system.stages().iter().filter(|s| s.id >= 0).count();
let models = (0..n_hydros)
.map(|h| {
if evap_indices.contains(&h) {
let coefficients = (0..n_stages)
.map(|s| LinearizedEvaporation {
k_evap0: k_evap0_per_stage
.get(s)
.copied()
.unwrap_or(k_evap0_per_stage.first().copied().unwrap_or(0.0)),
k_evap_v: 0.0,
})
.collect();
EvaporationModel::Linearized {
coefficients,
reference_volumes_hm3: vec![100.0; n_stages],
}
} else {
EvaporationModel::None
}
})
.collect();
EvaporationModelSet::new(models)
}
#[test]
fn evap_zero_hydros_layout_unchanged() {
let system = one_hydro_system(1, 0);
let no_evap = default_evaporation(&system);
let with_evap = build_stage_templates(
&system,
&no_penalty_config(),
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&default_production(&system),
&no_evap,
)
.expect("no evaporation ok");
let baseline = build_stage_templates(
&system,
&no_penalty_config(),
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&default_production(&system),
&EvaporationModelSet::new(vec![EvaporationModel::None]),
)
.expect("none evaporation ok");
assert_eq!(
with_evap.templates[0].num_cols, baseline.templates[0].num_cols,
"num_cols must match with zero evaporation hydros"
);
assert_eq!(
with_evap.templates[0].num_rows, baseline.templates[0].num_rows,
"num_rows must match with zero evaporation hydros"
);
}
#[test]
fn evap_two_hydros_increases_cols_and_rows() {
let system1 = one_hydro_system(1, 0);
let baseline = build_stage_templates(
&system1,
&no_penalty_config(),
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&default_production(&system1),
&EvaporationModelSet::new(vec![EvaporationModel::None]),
)
.expect("no evaporation baseline ok");
let evap = evap_set_for_system(&system1, &[0], &[1.5]);
let with_evap = build_stage_templates(
&system1,
&no_penalty_config(),
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&default_production(&system1),
&evap,
)
.expect("1 evaporation hydro ok");
let base_cols = baseline.templates[0].num_cols;
let base_rows = baseline.templates[0].num_rows;
let evap_cols = with_evap.templates[0].num_cols;
let evap_rows = with_evap.templates[0].num_rows;
assert_eq!(
evap_cols,
base_cols + 3,
"1 evap hydro must add exactly 3 columns (Q_ev, f_evap_plus, f_evap_minus)"
);
assert_eq!(
evap_rows,
base_rows + 1,
"1 evap hydro must add exactly 1 row (evaporation equality constraint)"
);
}
#[test]
fn evap_row_bounds_equality_at_k_evap0() {
let system = one_hydro_system(1, 0);
let k_evap0 = 1.5_f64;
let evap = evap_set_for_system(&system, &[0], &[k_evap0]);
let result = build_stage_templates(
&system,
&no_penalty_config(),
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&default_production(&system),
&evap,
)
.expect("evaporation system ok");
let t = &result.templates[0];
let evap_row = t.num_rows - 1;
assert_eq!(
t.row_lower[evap_row], k_evap0,
"evaporation row_lower must equal k_evap0 = {k_evap0}, got {}",
t.row_lower[evap_row]
);
assert_eq!(
t.row_upper[evap_row], k_evap0,
"evaporation row_upper must equal k_evap0 = {k_evap0}, got {}",
t.row_upper[evap_row]
);
}
#[test]
fn evap_col_bounds_and_objective() {
let system = one_hydro_system(1, 0);
let evap = evap_set_for_system(&system, &[0], &[1.5]);
let result = build_stage_templates(
&system,
&no_penalty_config(),
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&default_production(&system),
&evap,
)
.expect("evaporation system ok");
let t = &result.templates[0];
let col_q_ev = t.num_cols - 4;
let col_f_plus = t.num_cols - 3;
let col_f_minus = t.num_cols - 2;
for &col in &[col_q_ev, col_f_plus, col_f_minus] {
assert_eq!(
t.col_lower[col], 0.0,
"evap column {col} lower bound must be 0.0, got {}",
t.col_lower[col]
);
assert_eq!(
t.objective[col], 0.0,
"evap column {col} objective must be 0.0 (ticket-013 sets violation cost), got {}",
t.objective[col]
);
}
let expected_q_ev_bound = 1.5 * Q_EV_SAFETY_MARGIN;
assert!(
(t.col_upper[col_q_ev] - expected_q_ev_bound).abs() < 1e-12,
"Q_ev upper bound must be {expected_q_ev_bound}, got {}",
t.col_upper[col_q_ev]
);
for &col in &[col_f_plus, col_f_minus] {
assert!(
t.col_upper[col].is_infinite() && t.col_upper[col] > 0.0,
"evap slack column {col} upper bound must be +inf, got {}",
t.col_upper[col]
);
}
}
fn evap_set_with_k_evap_v(
system: &cobre_core::System,
evap_indices: &[usize],
k_evap0: f64,
k_evap_v: f64,
) -> EvaporationModelSet {
let n_hydros = system.hydros().len();
let n_stages = system.stages().iter().filter(|s| s.id >= 0).count();
let models = (0..n_hydros)
.map(|h| {
if evap_indices.contains(&h) {
let coefficients = (0..n_stages)
.map(|_| LinearizedEvaporation { k_evap0, k_evap_v })
.collect();
EvaporationModel::Linearized {
coefficients,
reference_volumes_hm3: vec![100.0; n_stages],
}
} else {
EvaporationModel::None
}
})
.collect();
EvaporationModelSet::new(models)
}
#[allow(clippy::cast_sign_loss)] fn entries_for_col(t: &StageTemplate, col: usize) -> Vec<(usize, f64)> {
let start = t.col_starts[col] as usize;
let end = t.col_starts[col + 1] as usize;
(start..end)
.map(|i| (t.row_indices[i] as usize, t.values[i]))
.collect()
}
#[test]
fn evap_csc_entries_one_hydro_correct_coefficients() {
let system = one_hydro_system(1, 0);
let k_evap_v = 0.02_f64;
let evap = evap_set_with_k_evap_v(&system, &[0], 1.5, k_evap_v);
let result = build_stage_templates(
&system,
&no_penalty_config(),
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&default_production(&system),
&evap,
)
.expect("evaporation system ok");
let t = &result.templates[0];
let col_q_ev = t.num_cols - 4;
let col_f_plus = t.num_cols - 3;
let col_f_minus = t.num_cols - 2;
let evap_row = t.num_rows - 1;
let water_balance_row = 2_usize;
let zeta = 744.0 * (3_600.0 / 1_000_000.0);
let entries_q_ev = entries_for_col(t, col_q_ev);
assert_eq!(
entries_q_ev.len(),
2,
"Q_ev column must have exactly 2 entries (water balance + evap constraint), got {entries_q_ev:?}"
);
assert_eq!(
entries_q_ev[0].0, water_balance_row,
"Q_ev first entry must be at water balance row"
);
assert!(
(entries_q_ev[0].1 - zeta).abs() < 1e-12,
"Q_ev water balance coefficient must be +zeta={zeta}, got {}",
entries_q_ev[0].1
);
assert_eq!(
entries_q_ev[1].0, evap_row,
"Q_ev second entry must be at evap_row"
);
assert!(
(entries_q_ev[1].1 - 1.0).abs() < 1e-12,
"Q_ev evap constraint coefficient must be +1.0, got {}",
entries_q_ev[1].1
);
let entries_f_plus = entries_for_col(t, col_f_plus);
assert_eq!(
entries_f_plus.len(),
1,
"f_plus column must have exactly 1 entry, got {entries_f_plus:?}"
);
assert_eq!(
entries_f_plus[0].0, evap_row,
"f_plus entry must be at evap_row"
);
assert!(
(entries_f_plus[0].1 - 1.0).abs() < 1e-12,
"f_plus coefficient must be +1.0, got {}",
entries_f_plus[0].1
);
let entries_f_minus = entries_for_col(t, col_f_minus);
assert_eq!(
entries_f_minus.len(),
1,
"f_minus column must have exactly 1 entry, got {entries_f_minus:?}"
);
assert_eq!(
entries_f_minus[0].0, evap_row,
"f_minus entry must be at evap_row"
);
assert!(
(entries_f_minus[0].1 - (-1.0)).abs() < 1e-12,
"f_minus coefficient must be -1.0, got {}",
entries_f_minus[0].1
);
let expected_coeff = -k_evap_v / 2.0;
let entry_v = entries_for_col(t, 0)
.into_iter()
.find(|&(r, _)| r == evap_row)
.expect("v column must have an entry at evap_row");
assert!(
(entry_v.1 - expected_coeff).abs() < 1e-12,
"v coefficient must be {expected_coeff}, got {}",
entry_v.1
);
let col_v_in = 2;
let entry_v_in = entries_for_col(t, col_v_in)
.into_iter()
.find(|&(r, _)| r == evap_row)
.expect("v_in column must have an entry at evap_row");
assert!(
(entry_v_in.1 - expected_coeff).abs() < 1e-12,
"v_in coefficient must be {expected_coeff}, got {}",
entry_v_in.1
);
}
#[test]
fn evap_csc_entries_coefficient_scaling() {
let system = one_hydro_system(1, 0);
let k_evap_v = 0.04_f64;
let evap = evap_set_with_k_evap_v(&system, &[0], 0.0, k_evap_v);
let result = build_stage_templates(
&system,
&no_penalty_config(),
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&default_production(&system),
&evap,
)
.expect("evaporation system ok");
let t = &result.templates[0];
let evap_row = t.num_rows - 1;
let expected_coeff = -k_evap_v / 2.0;
let entry_v = entries_for_col(t, 0)
.into_iter()
.find(|&(r, _)| r == evap_row)
.expect("v column must have evap_row entry");
assert!(
(entry_v.1 - expected_coeff).abs() < 1e-12,
"v coefficient: expected {expected_coeff}, got {}",
entry_v.1
);
let col_v_in = 2;
let entry_v_in = entries_for_col(t, col_v_in)
.into_iter()
.find(|&(r, _)| r == evap_row)
.expect("v_in column must have evap_row entry");
assert!(
(entry_v_in.1 - expected_coeff).abs() < 1e-12,
"v_in coefficient: expected {expected_coeff}, got {}",
entry_v_in.1
);
}
#[test]
fn evap_csc_entries_zero_hydros_no_op() {
let system = one_hydro_system(1, 0);
let no_evap = default_evaporation(&system);
let result = build_stage_templates(
&system,
&no_penalty_config(),
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&default_production(&system),
&no_evap,
)
.expect("no evaporation ok");
let baseline = build_stage_templates(
&system,
&no_penalty_config(),
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&default_production(&system),
&EvaporationModelSet::new(vec![EvaporationModel::None]),
)
.expect("none evaporation ok");
assert_eq!(
result.templates[0].num_nz, baseline.templates[0].num_nz,
"num_nz must be identical with zero evaporation hydros"
);
}
#[test]
fn evap_csc_entries_two_hydros_independent_rows() {
let (system, production) = four_hydro_mixed_system();
let n_stages = system.stages().iter().filter(|s| s.id >= 0).count();
let models = vec![
EvaporationModel::Linearized {
coefficients: vec![
LinearizedEvaporation {
k_evap0: 1.0,
k_evap_v: 0.02,
};
n_stages
],
reference_volumes_hm3: vec![100.0; n_stages],
},
EvaporationModel::Linearized {
coefficients: vec![
LinearizedEvaporation {
k_evap0: 2.0,
k_evap_v: 0.06,
};
n_stages
],
reference_volumes_hm3: vec![100.0; n_stages],
},
EvaporationModel::None,
EvaporationModel::None,
];
let evap = EvaporationModelSet::new(models);
let result = build_stage_templates(
&system,
&no_penalty_config(),
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&production,
&evap,
)
.expect("2-evap-hydro system ok");
let t = &result.templates[0];
let evap_row_0 = t.num_rows - 2;
let evap_row_1 = t.num_rows - 1;
let entry_v_h0 = entries_for_col(t, 0)
.into_iter()
.find(|&(r, _)| r == evap_row_0)
.expect("hydro 0 v col entry");
assert!(
(entry_v_h0.1 - (-0.01)).abs() < 1e-12,
"hydro 0 v: expected -0.01, got {}",
entry_v_h0.1
);
let entry_v_h1 = entries_for_col(t, 1)
.into_iter()
.find(|&(r, _)| r == evap_row_1)
.expect("hydro 1 v col entry");
assert!(
(entry_v_h1.1 - (-0.03)).abs() < 1e-12,
"hydro 1 v: expected -0.03, got {}",
entry_v_h1.1
);
assert!((t.row_lower[evap_row_0] - 1.0).abs() < 1e-12);
assert!((t.row_lower[evap_row_1] - 2.0).abs() < 1e-12);
}
#[test]
fn evap_csc_entries_zero_k_evap_v_produces_zero_volume_coefficients() {
let system = one_hydro_system(1, 0);
let evap = evap_set_with_k_evap_v(&system, &[0], 2.0, 0.0);
let result = build_stage_templates(
&system,
&no_penalty_config(),
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&default_production(&system),
&evap,
)
.expect("evaporation system ok");
let t = &result.templates[0];
let evap_row = t.num_rows - 1;
let entry_v = entries_for_col(t, 0)
.into_iter()
.find(|&(r, _)| r == evap_row)
.expect("v column must have evap_row entry");
assert!(
entry_v.1.abs() < 1e-12,
"v coefficient must be 0.0 when k_evap_v=0, got {}",
entry_v.1
);
let col_v_in = 2;
let entry_v_in = entries_for_col(t, col_v_in)
.into_iter()
.find(|&(r, _)| r == evap_row)
.expect("v_in column must have evap_row entry");
assert!(
entry_v_in.1.abs() < 1e-12,
"v_in coefficient must be 0.0 when k_evap_v=0, got {}",
entry_v_in.1
);
}
#[test]
#[allow(clippy::cast_sign_loss)]
fn evap_water_balance_one_hydro_coefficient_is_zeta() {
let system = one_hydro_system(1, 0);
let evap = evap_set_with_k_evap_v(&system, &[0], 0.0, 0.0);
let result = build_stage_templates(
&system,
&no_penalty_config(),
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&default_production(&system),
&evap,
)
.expect("evaporation system ok");
let t = &result.templates[0];
let water_balance_row = 2_usize;
let col_q_ev = t.num_cols - 4;
let entries = entries_for_col(t, col_q_ev);
let entry = entries
.iter()
.find(|&&(r, _)| r == water_balance_row)
.copied()
.expect("Q_ev column must have an entry in the water balance row");
let zeta = 744.0_f64 * (3_600.0 / 1_000_000.0);
assert!(
(entry.1 - zeta).abs() < 1e-12,
"Q_ev water balance coefficient must be +zeta={zeta}, got {}",
entry.1
);
}
#[test]
#[allow(clippy::cast_sign_loss, clippy::too_many_lines)]
fn evap_water_balance_only_second_hydro_has_evap() {
use chrono::NaiveDate;
use cobre_core::entities::hydro::{Hydro, HydroGenerationModel, HydroPenalties};
use cobre_core::scenario::{InflowModel, LoadModel};
use cobre_core::temporal::{
Block, BlockMode, NoiseMethod, ScenarioSourceConfig, Stage as CStage, StageRiskConfig,
StageStateConfig,
};
let bus = Bus {
id: EntityId(1),
name: "B1".to_string(),
deficit_segments: vec![DeficitSegment {
depth_mw: None,
cost_per_mwh: 500.0,
}],
excess_cost: 0.0,
};
let hp = HydroPenalties {
spillage_cost: 0.01,
diversion_cost: 0.0,
fpha_turbined_cost: 0.0,
storage_violation_below_cost: 0.0,
filling_target_violation_cost: 0.0,
turbined_violation_below_cost: 0.0,
outflow_violation_below_cost: 0.0,
outflow_violation_above_cost: 0.0,
generation_violation_below_cost: 0.0,
evaporation_violation_cost: 0.0,
water_withdrawal_violation_cost: 0.0,
};
let make_h = |id: i32| Hydro {
id: EntityId(id),
name: format!("H{id}"),
bus_id: EntityId(1),
downstream_id: None,
entry_stage_id: None,
exit_stage_id: None,
min_storage_hm3: 0.0,
max_storage_hm3: 200.0,
min_outflow_m3s: 0.0,
max_outflow_m3s: None,
generation_model: HydroGenerationModel::ConstantProductivity {
productivity_mw_per_m3s: 2.5,
},
min_turbined_m3s: 0.0,
max_turbined_m3s: 100.0,
min_generation_mw: 0.0,
max_generation_mw: 250.0,
tailrace: None,
hydraulic_losses: None,
efficiency: None,
evaporation_coefficients_mm: None,
evaporation_reference_volumes_hm3: None,
diversion: None,
filling: None,
penalties: hp,
};
let hydros = vec![make_h(2), make_h(3)];
let stages = vec![CStage {
index: 0,
id: 0,
start_date: NaiveDate::from_ymd_opt(2024, 1, 1).unwrap(),
end_date: NaiveDate::from_ymd_opt(2024, 2, 1).unwrap(),
season_id: None,
blocks: vec![Block {
index: 0,
name: "S".to_string(),
duration_hours: 744.0,
}],
block_mode: BlockMode::Parallel,
state_config: StageStateConfig {
storage: true,
inflow_lags: false,
},
risk_config: StageRiskConfig::Expectation,
scenario_config: ScenarioSourceConfig {
branching_factor: 1,
noise_method: NoiseMethod::Saa,
},
}];
let inflow_models: Vec<InflowModel> = hydros
.iter()
.map(|h| InflowModel {
hydro_id: h.id,
stage_id: 0,
mean_m3s: 80.0,
std_m3s: 20.0,
ar_coefficients: vec![],
residual_std_ratio: 1.0,
})
.collect();
let load_models = vec![LoadModel {
bus_id: EntityId(1),
stage_id: 0,
mean_mw: 200.0,
std_mw: 0.0,
}];
let bounds = ResolvedBounds::new(
2,
0,
0,
0,
0,
1,
default_hydro_bounds(),
ThermalStageBounds {
min_generation_mw: 0.0,
max_generation_mw: 0.0,
},
LineStageBounds {
direct_mw: 0.0,
reverse_mw: 0.0,
},
PumpingStageBounds {
min_flow_m3s: 0.0,
max_flow_m3s: 0.0,
},
ContractStageBounds {
min_mw: 0.0,
max_mw: 0.0,
price_per_mwh: 0.0,
},
);
let penalties = ResolvedPenalties::new(
2,
1,
0,
0,
1,
default_hydro_penalties(),
BusStagePenalties { excess_cost: 0.0 },
LineStagePenalties { exchange_cost: 0.0 },
NcsStagePenalties {
curtailment_cost: 0.0,
},
);
let system = SystemBuilder::new()
.buses(vec![bus])
.hydros(hydros)
.stages(stages)
.inflow_models(inflow_models)
.load_models(load_models)
.bounds(bounds)
.penalties(penalties)
.build()
.expect("2-hydro system ok");
let evap = evap_set_with_k_evap_v(&system, &[1], 0.0, 0.0);
let result = build_stage_templates(
&system,
&no_penalty_config(),
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&default_production(&system),
&evap,
)
.expect("2-hydro evap system ok");
let t = &result.templates[0];
let water_balance_row_h0 = 4_usize;
let water_balance_row_h1 = 5_usize;
let col_q_ev_h1 = t.num_cols - 5;
let entries_h1 = entries_for_col(t, col_q_ev_h1);
let found_h1 = entries_h1
.iter()
.find(|&&(r, _)| r == water_balance_row_h1)
.copied();
assert!(
found_h1.is_some(),
"Q_ev for hydro 1 must have an entry in water balance row {water_balance_row_h1}"
);
let zeta = 744.0_f64 * (3_600.0 / 1_000_000.0);
assert!(
(found_h1.unwrap().1 - zeta).abs() < 1e-12,
"Q_ev (h1) water balance coefficient must be +zeta={zeta}, got {}",
found_h1.unwrap().1
);
let found_h0 = entries_h1.iter().any(|&(r, _)| r == water_balance_row_h0);
assert!(
!found_h0,
"Q_ev for hydro 1 must not appear in hydro 0's water balance row"
);
}
#[test]
fn evap_water_balance_zero_hydros_no_op() {
let system = one_hydro_system(1, 0);
let no_evap = EvaporationModelSet::new(vec![EvaporationModel::None]);
let result = build_stage_templates(
&system,
&no_penalty_config(),
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&default_production(&system),
&no_evap,
)
.expect("no evaporation ok");
let baseline = build_stage_templates(
&system,
&no_penalty_config(),
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&default_production(&system),
&default_evaporation(&system),
)
.expect("default evaporation ok");
assert_eq!(
result.templates[0].num_nz, baseline.templates[0].num_nz,
"num_nz must be identical with zero evaporation hydros (no water balance entries added)"
);
}
#[allow(
clippy::cast_possible_truncation,
clippy::cast_possible_wrap,
clippy::too_many_lines
)]
fn evap_hydro_system_with_violation_cost(
block_hours: f64,
evaporation_violation_cost: f64,
) -> cobre_core::System {
use chrono::NaiveDate;
use cobre_core::entities::hydro::{Hydro, HydroGenerationModel, HydroPenalties};
use cobre_core::scenario::{InflowModel, LoadModel};
use cobre_core::temporal::{
Block, BlockMode, NoiseMethod, ScenarioSourceConfig, Stage, StageRiskConfig,
StageStateConfig,
};
let bus = Bus {
id: EntityId(1),
name: "B1".to_string(),
deficit_segments: vec![DeficitSegment {
depth_mw: None,
cost_per_mwh: 500.0,
}],
excess_cost: 0.0,
};
let hydro = Hydro {
id: EntityId(2),
name: "H1".to_string(),
bus_id: EntityId(1),
downstream_id: None,
entry_stage_id: None,
exit_stage_id: None,
min_storage_hm3: 0.0,
max_storage_hm3: 2_000.0,
min_outflow_m3s: 0.0,
max_outflow_m3s: None,
generation_model: HydroGenerationModel::ConstantProductivity {
productivity_mw_per_m3s: 2.5,
},
min_turbined_m3s: 0.0,
max_turbined_m3s: 100.0,
min_generation_mw: 0.0,
max_generation_mw: 250.0,
tailrace: None,
hydraulic_losses: None,
efficiency: None,
evaporation_coefficients_mm: None,
evaporation_reference_volumes_hm3: None,
diversion: None,
filling: None,
penalties: HydroPenalties {
spillage_cost: 0.01,
diversion_cost: 0.0,
fpha_turbined_cost: 0.0,
storage_violation_below_cost: 0.0,
filling_target_violation_cost: 0.0,
turbined_violation_below_cost: 0.0,
outflow_violation_below_cost: 0.0,
outflow_violation_above_cost: 0.0,
generation_violation_below_cost: 0.0,
evaporation_violation_cost,
water_withdrawal_violation_cost: 0.0,
},
};
let stages = vec![Stage {
index: 0,
id: 0,
start_date: NaiveDate::from_ymd_opt(2024, 1, 1).unwrap(),
end_date: NaiveDate::from_ymd_opt(2024, 2, 1).unwrap(),
season_id: None,
blocks: vec![Block {
index: 0,
name: "S".to_string(),
duration_hours: block_hours,
}],
block_mode: BlockMode::Parallel,
state_config: StageStateConfig {
storage: true,
inflow_lags: false,
},
risk_config: StageRiskConfig::Expectation,
scenario_config: ScenarioSourceConfig {
branching_factor: 1,
noise_method: NoiseMethod::Saa,
},
}];
let inflow_models = vec![InflowModel {
hydro_id: EntityId(2),
stage_id: 0,
mean_m3s: 50.0,
std_m3s: 10.0,
ar_coefficients: vec![],
residual_std_ratio: 1.0,
}];
let load_models = vec![LoadModel {
bus_id: EntityId(1),
stage_id: 0,
mean_mw: 100.0,
std_mw: 0.0,
}];
let bounds = ResolvedBounds::new(
1,
0,
0,
0,
0,
1,
HydroStageBounds {
min_storage_hm3: 0.0,
max_storage_hm3: 2_000.0,
min_turbined_m3s: 0.0,
max_turbined_m3s: 100.0,
min_outflow_m3s: 0.0,
max_outflow_m3s: None,
min_generation_mw: 0.0,
max_generation_mw: 250.0,
max_diversion_m3s: None,
filling_inflow_m3s: 0.0,
water_withdrawal_m3s: 0.0,
},
ThermalStageBounds {
min_generation_mw: 0.0,
max_generation_mw: 0.0,
},
LineStageBounds {
direct_mw: 0.0,
reverse_mw: 0.0,
},
PumpingStageBounds {
min_flow_m3s: 0.0,
max_flow_m3s: 0.0,
},
ContractStageBounds {
min_mw: 0.0,
max_mw: 0.0,
price_per_mwh: 0.0,
},
);
let penalties = ResolvedPenalties::new(
1,
1,
0,
0,
1,
HydroStagePenalties {
spillage_cost: 0.01,
diversion_cost: 0.0,
fpha_turbined_cost: 0.0,
storage_violation_below_cost: 0.0,
filling_target_violation_cost: 0.0,
turbined_violation_below_cost: 0.0,
outflow_violation_below_cost: 0.0,
outflow_violation_above_cost: 0.0,
generation_violation_below_cost: 0.0,
evaporation_violation_cost,
water_withdrawal_violation_cost: 0.0,
},
BusStagePenalties { excess_cost: 0.0 },
LineStagePenalties { exchange_cost: 0.0 },
NcsStagePenalties {
curtailment_cost: 0.0,
},
);
SystemBuilder::new()
.buses(vec![bus])
.hydros(vec![hydro])
.stages(stages)
.inflow_models(inflow_models)
.load_models(load_models)
.bounds(bounds)
.penalties(penalties)
.build()
.expect("evap_hydro_system_with_violation_cost: valid")
}
#[test]
fn evap_violation_cost_applied_to_slack_columns() {
let system = evap_hydro_system_with_violation_cost(730.0, 500.0);
let evap = evap_set_with_k_evap_v(&system, &[0], 1.0, 0.02);
let result = build_stage_templates(
&system,
&no_penalty_config(),
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&default_production(&system),
&evap,
)
.expect("evap violation cost system builds ok");
let t = &result.templates[0];
let col_q_ev = t.num_cols - 4;
let col_f_plus = t.num_cols - 3;
let col_f_minus = t.num_cols - 2;
let expected_base = 500.0 * 730.0 / COST_SCALE_FACTOR;
assert!(
t.objective[col_q_ev].abs() < 1e-12,
"Q_ev column objective must be 0.0 (evaporation flow itself has no cost), got {}",
t.objective[col_q_ev]
);
assert!(
(t.objective[col_f_plus] - expected_base).abs() < 1e-12,
"f_evap_plus objective: expected {expected_base}, got {}",
t.objective[col_f_plus]
);
let expected_minus = expected_base * OVER_EVAPORATION_COST_MULTIPLIER;
assert!(
(t.objective[col_f_minus] - expected_minus).abs() < 1e-6,
"f_evap_minus objective: expected {expected_minus} (100x f_plus), got {}",
t.objective[col_f_minus]
);
}
#[test]
fn evap_q_ev_objective_is_zero() {
let system = evap_hydro_system_with_violation_cost(730.0, 500.0);
let evap = evap_set_with_k_evap_v(&system, &[0], 0.0, 0.0);
let result = build_stage_templates(
&system,
&no_penalty_config(),
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&default_production(&system),
&evap,
)
.expect("evap system with zero k_evap builds ok");
let t = &result.templates[0];
let col_q_ev = t.num_cols - 4;
assert!(
t.objective[col_q_ev].abs() < 1e-12,
"Q_ev objective must be 0.0, got {}",
t.objective[col_q_ev]
);
}
#[test]
fn evap_lp_solvable_and_q_ev_nonnegative() {
use cobre_solver::{HighsSolver, RowBatch, SolverInterface};
let system = evap_hydro_system_with_violation_cost(730.0, 500.0);
let evap = evap_set_with_k_evap_v(&system, &[0], 1.0, 0.02);
let result = build_stage_templates(
&system,
&no_penalty_config(),
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&default_production(&system),
&evap,
)
.expect("evap system template build must succeed");
let template = &result.templates[0];
let mut solver = HighsSolver::new().expect("HighsSolver::new must succeed");
solver.load_model(template);
let empty_cuts = RowBatch {
num_rows: 0,
row_starts: vec![0_i32],
col_indices: vec![],
values: vec![],
row_lower: vec![],
row_upper: vec![],
};
solver.add_rows(&empty_cuts);
let v_in = 1_000.0_f64;
solver.set_row_bounds(&[0], &[v_in], &[v_in]);
let view = solver
.solve()
.expect("evaporation LP must be feasible and optimal");
let col_q_ev = template.num_cols - 4;
let q_ev = view.primal[col_q_ev];
assert!(
q_ev >= -1e-8,
"Q_ev must be non-negative after solving, got {q_ev}"
);
}
#[test]
fn evap_violation_slacks_near_zero_feasible_constraint() {
use cobre_solver::{HighsSolver, RowBatch, SolverInterface};
let system = evap_hydro_system_with_violation_cost(730.0, 500.0);
let evap = evap_set_with_k_evap_v(&system, &[0], 1.0, 0.02);
let result = build_stage_templates(
&system,
&no_penalty_config(),
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&default_production(&system),
&evap,
)
.expect("evap system template build must succeed");
let template = &result.templates[0];
let mut solver = HighsSolver::new().expect("HighsSolver::new must succeed");
solver.load_model(template);
let empty_cuts = RowBatch {
num_rows: 0,
row_starts: vec![0_i32],
col_indices: vec![],
values: vec![],
row_lower: vec![],
row_upper: vec![],
};
solver.add_rows(&empty_cuts);
let v_in = 1_000.0_f64;
solver.set_row_bounds(&[0], &[v_in], &[v_in]);
let view = solver
.solve()
.expect("evaporation LP must be feasible and optimal");
let col_f_plus = template.num_cols - 3;
let col_f_minus = template.num_cols - 2;
let f_plus = view.primal[col_f_plus];
let f_minus = view.primal[col_f_minus];
assert!(
f_plus.abs() < 1e-6,
"f_evap_plus slack must be near zero (constraint satisfied without violation), got {f_plus}"
);
assert!(
f_minus.abs() < 1e-6,
"f_evap_minus slack must be near zero (constraint satisfied without violation), got {f_minus}"
);
}
#[test]
fn evap_storage_fixing_dual_differs_from_no_evaporation() {
use cobre_solver::{HighsSolver, RowBatch, SolverInterface};
let system_evap = evap_hydro_system_with_violation_cost(730.0, 500.0);
let evap = evap_set_with_k_evap_v(&system_evap, &[0], 1.0, 0.02);
let evap_result = build_stage_templates(
&system_evap,
&no_penalty_config(),
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&default_production(&system_evap),
&evap,
)
.expect("evap system template build must succeed");
let system_base = one_hydro_system(1, 0);
let no_evap = EvaporationModelSet::new(vec![EvaporationModel::None]);
let base_result = build_stage_templates(
&system_base,
&no_penalty_config(),
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&default_production(&system_base),
&no_evap,
)
.expect("baseline system template build must succeed");
let solve_and_get_storage_dual = |template: &cobre_solver::StageTemplate| -> f64 {
let mut solver = HighsSolver::new().expect("HighsSolver::new must succeed");
solver.load_model(template);
let empty_cuts = RowBatch {
num_rows: 0,
row_starts: vec![0_i32],
col_indices: vec![],
values: vec![],
row_lower: vec![],
row_upper: vec![],
};
solver.add_rows(&empty_cuts);
let v_in = 1_000.0_f64;
solver.set_row_bounds(&[0], &[v_in], &[v_in]);
let view = solver.solve().expect("LP must solve to optimal");
view.dual[0]
};
let evap_dual = solve_and_get_storage_dual(&evap_result.templates[0]);
let base_dual = solve_and_get_storage_dual(&base_result.templates[0]);
assert_ne!(
(evap_dual * 1e6).round(),
(base_dual * 1e6).round(),
"storage-fixing dual must differ between evaporation ({evap_dual}) and \
no-evaporation ({base_dual}) configurations"
);
}
#[test]
fn evap_bound_prevents_dump_valve() {
use cobre_solver::{HighsSolver, RowBatch, SolverInterface};
let system = evap_hydro_system_with_violation_cost(730.0, 500.0);
let evap = evap_set_with_k_evap_v(&system, &[0], 2.0, 0.0001);
let result = build_stage_templates(
&system,
&no_penalty_config(),
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&default_production(&system),
&evap,
)
.expect("evap dump valve test: template build must succeed");
let template = &result.templates[0];
let mut solver = HighsSolver::new().expect("HighsSolver::new must succeed");
solver.load_model(template);
let empty_cuts = RowBatch {
num_rows: 0,
row_starts: vec![0_i32],
col_indices: vec![],
values: vec![],
row_lower: vec![],
row_upper: vec![],
};
solver.add_rows(&empty_cuts);
let v_in = 2_000.0_f64;
solver.set_row_bounds(&[0], &[v_in], &[v_in]);
let zeta = 730.0 * 3600.0 / 1e6;
let high_inflow_rhs = zeta * 500.0;
solver.set_row_bounds(&[2], &[high_inflow_rhs], &[high_inflow_rhs]);
let view = solver
.solve()
.expect("evap dump valve LP must be feasible and optimal");
let col_spillage = 5;
let col_q_ev = template.num_cols - 4;
let col_f_minus = template.num_cols - 2;
let q_ev = view.primal[col_q_ev];
let f_minus = view.primal[col_f_minus];
let spillage = view.primal[col_spillage];
let q_ev_max = (2.0 + 0.0001 * 2_000.0) * Q_EV_SAFETY_MARGIN;
assert!(
q_ev <= q_ev_max + 1e-8,
"Q_ev must be bounded by physical limit {q_ev_max}, got {q_ev}"
);
assert!(
f_minus < 1e-6,
"f_minus (over-evaporation) must be near zero, got {f_minus}"
);
assert!(
spillage > 1e-6,
"spillage must be positive when excess water needs dumping, got {spillage}"
);
}
#[allow(clippy::cast_possible_truncation, clippy::cast_possible_wrap)]
fn multi_segment_system(buses: Vec<Bus>, block_hours: f64) -> cobre_core::System {
use chrono::NaiveDate;
use cobre_core::scenario::LoadModel;
use cobre_core::temporal::{
Block, BlockMode, ScenarioSourceConfig, Stage, StageRiskConfig, StageStateConfig,
};
let n_buses = buses.len();
let load_models: Vec<LoadModel> = buses
.iter()
.map(|b| LoadModel {
bus_id: b.id,
stage_id: 0,
mean_mw: 0.0,
std_mw: 0.0,
})
.collect();
let stage = Stage {
index: 0,
id: 0,
start_date: NaiveDate::from_ymd_opt(2024, 1, 1).unwrap(),
end_date: NaiveDate::from_ymd_opt(2024, 2, 1).unwrap(),
season_id: None,
blocks: vec![Block {
index: 0,
name: "S".to_string(),
duration_hours: block_hours,
}],
block_mode: BlockMode::Parallel,
state_config: StageStateConfig {
storage: false,
inflow_lags: false,
},
risk_config: StageRiskConfig::Expectation,
scenario_config: ScenarioSourceConfig {
branching_factor: 1,
noise_method: cobre_core::temporal::NoiseMethod::Saa,
},
};
let bounds = ResolvedBounds::new(
0,
0,
0,
0,
n_buses,
1,
default_hydro_bounds(),
ThermalStageBounds {
min_generation_mw: 0.0,
max_generation_mw: 0.0,
},
LineStageBounds {
direct_mw: 0.0,
reverse_mw: 0.0,
},
PumpingStageBounds {
min_flow_m3s: 0.0,
max_flow_m3s: 0.0,
},
ContractStageBounds {
min_mw: 0.0,
max_mw: 0.0,
price_per_mwh: 0.0,
},
);
let penalties = ResolvedPenalties::new(
0,
n_buses,
0,
0,
1,
default_hydro_penalties(),
BusStagePenalties { excess_cost: 0.0 },
LineStagePenalties { exchange_cost: 0.0 },
NcsStagePenalties {
curtailment_cost: 0.0,
},
);
SystemBuilder::new()
.buses(buses)
.stages(vec![stage])
.load_models(load_models)
.bounds(bounds)
.penalties(penalties)
.build()
.expect("multi_segment_system: valid")
}
#[test]
#[allow(clippy::too_many_lines)]
fn test_multi_segment_deficit_column_count() {
use chrono::NaiveDate;
use cobre_core::scenario::LoadModel;
use cobre_core::temporal::{
Block, BlockMode, ScenarioSourceConfig, Stage, StageRiskConfig, StageStateConfig,
};
let bus0 = Bus {
id: EntityId(1),
name: "Bus0".to_string(),
deficit_segments: vec![
DeficitSegment {
depth_mw: Some(10.0),
cost_per_mwh: 100.0,
},
DeficitSegment {
depth_mw: Some(20.0),
cost_per_mwh: 200.0,
},
DeficitSegment {
depth_mw: None,
cost_per_mwh: 5000.0,
},
],
excess_cost: 0.0,
};
let bus1 = Bus {
id: EntityId(2),
name: "Bus1".to_string(),
deficit_segments: vec![DeficitSegment {
depth_mw: None,
cost_per_mwh: 1000.0,
}],
excess_cost: 0.0,
};
let load_models = vec![
LoadModel {
bus_id: EntityId(1),
stage_id: 0,
mean_mw: 0.0,
std_mw: 0.0,
},
LoadModel {
bus_id: EntityId(2),
stage_id: 0,
mean_mw: 0.0,
std_mw: 0.0,
},
];
let stage = Stage {
index: 0,
id: 0,
start_date: NaiveDate::from_ymd_opt(2024, 1, 1).unwrap(),
end_date: NaiveDate::from_ymd_opt(2024, 2, 1).unwrap(),
season_id: None,
blocks: vec![
Block {
index: 0,
name: "B0".to_string(),
duration_hours: 360.0,
},
Block {
index: 1,
name: "B1".to_string(),
duration_hours: 360.0,
},
],
block_mode: BlockMode::Parallel,
state_config: StageStateConfig {
storage: false,
inflow_lags: false,
},
risk_config: StageRiskConfig::Expectation,
scenario_config: ScenarioSourceConfig {
branching_factor: 1,
noise_method: cobre_core::temporal::NoiseMethod::Saa,
},
};
let bounds = ResolvedBounds::new(
0,
0,
0,
0,
2,
1,
default_hydro_bounds(),
ThermalStageBounds {
min_generation_mw: 0.0,
max_generation_mw: 0.0,
},
LineStageBounds {
direct_mw: 0.0,
reverse_mw: 0.0,
},
PumpingStageBounds {
min_flow_m3s: 0.0,
max_flow_m3s: 0.0,
},
ContractStageBounds {
min_mw: 0.0,
max_mw: 0.0,
price_per_mwh: 0.0,
},
);
let penalties = ResolvedPenalties::new(
0,
2,
0,
0,
1,
default_hydro_penalties(),
BusStagePenalties { excess_cost: 0.0 },
LineStagePenalties { exchange_cost: 0.0 },
NcsStagePenalties {
curtailment_cost: 0.0,
},
);
let system = SystemBuilder::new()
.buses(vec![bus0, bus1])
.stages(vec![stage])
.load_models(load_models)
.bounds(bounds)
.penalties(penalties)
.build()
.expect("2-bus 2-block system: valid");
let result = build_stage_templates(
&system,
&no_penalty_config(),
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&default_production(&system),
&default_evaporation(&system),
)
.expect("build ok");
let t = &result.templates[0];
let col_deficit_start = 1_usize;
let max_deficit_segments = 3_usize;
let n_blks = 2_usize;
let n_buses = 2_usize;
let deficit_region = n_buses * max_deficit_segments * n_blks;
assert_eq!(
deficit_region, 12,
"deficit region must be B*S_max*K = 2*3*2 = 12"
);
let col_excess_start = col_deficit_start + deficit_region;
let excess_region = n_buses * n_blks; let expected_num_cols = col_excess_start + excess_region;
assert_eq!(
t.num_cols, expected_num_cols,
"num_cols must include expanded deficit region"
);
}
#[test]
fn test_multi_segment_deficit_bounds_and_objective() {
let bus = Bus {
id: EntityId(1),
name: "Bus0".to_string(),
deficit_segments: vec![
DeficitSegment {
depth_mw: Some(10.0),
cost_per_mwh: 500.0,
},
DeficitSegment {
depth_mw: None,
cost_per_mwh: 5000.0,
},
],
excess_cost: 0.0,
};
let block_hours = 730.0_f64;
let system = multi_segment_system(vec![bus], block_hours);
let result = build_stage_templates(
&system,
&no_penalty_config(),
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&default_production(&system),
&default_evaporation(&system),
)
.expect("build ok");
let t = &result.templates[0];
let col_seg0 = 1_usize;
let col_seg1 = 2_usize;
assert_eq!(
t.col_upper[col_seg0], 10.0,
"segment 0 upper bound must equal depth_mw = 10.0"
);
assert!(
t.col_upper[col_seg1].is_infinite() && t.col_upper[col_seg1] > 0.0,
"segment 1 upper bound must be +infinity (unbounded final segment)"
);
assert!(
(t.objective[col_seg0] - 500.0 * block_hours / COST_SCALE_FACTOR).abs() < 1e-12,
"segment 0 objective must be cost * block_hours / COST_SCALE_FACTOR = {} but got {}",
500.0 * block_hours / COST_SCALE_FACTOR,
t.objective[col_seg0]
);
assert!(
(t.objective[col_seg1] - 5000.0 * block_hours / COST_SCALE_FACTOR).abs() < 1e-12,
"segment 1 objective must be cost * block_hours / COST_SCALE_FACTOR = {} but got {}",
5000.0 * block_hours / COST_SCALE_FACTOR,
t.objective[col_seg1]
);
}
#[test]
fn test_single_segment_backward_compat() {
let cost = 1000.0_f64;
let block_hours = 744.0_f64;
let bus = Bus {
id: EntityId(1),
name: "Bus0".to_string(),
deficit_segments: vec![DeficitSegment {
depth_mw: None,
cost_per_mwh: cost,
}],
excess_cost: 0.0,
};
let system = multi_segment_system(vec![bus], block_hours);
let result = build_stage_templates(
&system,
&no_penalty_config(),
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&default_production(&system),
&default_evaporation(&system),
)
.expect("build ok");
let t = &result.templates[0];
let col_def = 1_usize;
assert!(
t.col_upper[col_def].is_infinite() && t.col_upper[col_def] > 0.0,
"single segment must be unbounded (None depth_mw)"
);
assert!(
(t.objective[col_def] - cost * block_hours / COST_SCALE_FACTOR).abs() < 1e-12,
"single-segment objective must be cost * block_hours / COST_SCALE_FACTOR"
);
let col_exc = 2_usize;
assert!(
t.col_upper[col_exc].is_infinite() && t.col_upper[col_exc] > 0.0,
"excess column must be unbounded"
);
}
#[test]
fn test_multi_segment_deficit_load_balance_coefficients() {
let bus = Bus {
id: EntityId(1),
name: "Bus0".to_string(),
deficit_segments: vec![
DeficitSegment {
depth_mw: Some(10.0),
cost_per_mwh: 500.0,
},
DeficitSegment {
depth_mw: None,
cost_per_mwh: 5000.0,
},
],
excess_cost: 0.0,
};
let system = multi_segment_system(vec![bus], 730.0);
let result = build_stage_templates(
&system,
&no_penalty_config(),
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&default_production(&system),
&default_evaporation(&system),
)
.expect("build ok");
let t = &result.templates[0];
let col_seg0 = 1_usize;
let col_seg1 = 2_usize;
let load_balance_row = 0_usize;
let entries_seg0 = entries_for_col(t, col_seg0);
assert_eq!(
entries_seg0.len(),
1,
"deficit segment 0 column must have exactly 1 CSC entry (load-balance row), got {entries_seg0:?}"
);
assert_eq!(
entries_seg0[0].0, load_balance_row,
"deficit segment 0 entry must be at the load-balance row {load_balance_row}, got row {}",
entries_seg0[0].0
);
assert!(
(entries_seg0[0].1 - 1.0).abs() < 1e-12,
"deficit segment 0 load-balance coefficient must be +1.0, got {}",
entries_seg0[0].1
);
let entries_seg1 = entries_for_col(t, col_seg1);
assert_eq!(
entries_seg1.len(),
1,
"deficit segment 1 column must have exactly 1 CSC entry (load-balance row), got {entries_seg1:?}"
);
assert_eq!(
entries_seg1[0].0, load_balance_row,
"deficit segment 1 entry must be at the load-balance row {load_balance_row}, got row {}",
entries_seg1[0].0
);
assert!(
(entries_seg1[0].1 - 1.0).abs() < 1e-12,
"deficit segment 1 load-balance coefficient must be +1.0, got {}",
entries_seg1[0].1
);
}
#[allow(
clippy::cast_possible_truncation,
clippy::cast_possible_wrap,
clippy::too_many_lines
)]
fn one_hydro_system_with_withdrawal(
n_stages: usize,
lag_order: usize,
water_withdrawal_m3s: f64,
water_withdrawal_violation_cost: f64,
) -> cobre_core::System {
use chrono::NaiveDate;
use cobre_core::entities::hydro::{Hydro, HydroGenerationModel, HydroPenalties};
use cobre_core::scenario::{InflowModel, LoadModel};
use cobre_core::temporal::{
Block, BlockMode, NoiseMethod, ScenarioSourceConfig, Stage, StageRiskConfig,
StageStateConfig,
};
let bus = Bus {
id: EntityId(1),
name: "B1".to_string(),
deficit_segments: vec![DeficitSegment {
depth_mw: None,
cost_per_mwh: 500.0,
}],
excess_cost: 0.0,
};
let hydro = Hydro {
id: EntityId(2),
name: "H1".to_string(),
bus_id: EntityId(1),
downstream_id: None,
entry_stage_id: None,
exit_stage_id: None,
min_storage_hm3: 0.0,
max_storage_hm3: 200.0,
min_outflow_m3s: 0.0,
max_outflow_m3s: None,
generation_model: HydroGenerationModel::ConstantProductivity {
productivity_mw_per_m3s: 2.5,
},
min_turbined_m3s: 0.0,
max_turbined_m3s: 100.0,
min_generation_mw: 0.0,
max_generation_mw: 250.0,
tailrace: None,
hydraulic_losses: None,
efficiency: None,
evaporation_coefficients_mm: None,
evaporation_reference_volumes_hm3: None,
diversion: None,
filling: None,
penalties: HydroPenalties {
spillage_cost: 0.01,
diversion_cost: 0.0,
fpha_turbined_cost: 0.0,
storage_violation_below_cost: 0.0,
filling_target_violation_cost: 0.0,
turbined_violation_below_cost: 0.0,
outflow_violation_below_cost: 0.0,
outflow_violation_above_cost: 0.0,
generation_violation_below_cost: 0.0,
evaporation_violation_cost: 0.0,
water_withdrawal_violation_cost,
},
};
let stages: Vec<Stage> = (0..n_stages)
.map(|i| Stage {
index: i,
id: i as i32,
start_date: NaiveDate::from_ymd_opt(2024, 1, 1).unwrap(),
end_date: NaiveDate::from_ymd_opt(2024, 2, 1).unwrap(),
season_id: None,
blocks: vec![Block {
index: 0,
name: "S".to_string(),
duration_hours: 744.0,
}],
block_mode: BlockMode::Parallel,
state_config: StageStateConfig {
storage: true,
inflow_lags: lag_order > 0,
},
risk_config: StageRiskConfig::Expectation,
scenario_config: ScenarioSourceConfig {
branching_factor: 1,
noise_method: NoiseMethod::Saa,
},
})
.collect();
let ar_coefficients: Vec<f64> = (0..lag_order).map(|_| 0.5).collect();
let inflow_models: Vec<InflowModel> = (0..n_stages)
.map(|i| InflowModel {
hydro_id: EntityId(2),
stage_id: i as i32,
mean_m3s: 80.0,
std_m3s: 20.0,
ar_coefficients: ar_coefficients.clone(),
residual_std_ratio: 1.0,
})
.collect();
let load_models: Vec<LoadModel> = (0..n_stages)
.map(|i| LoadModel {
bus_id: EntityId(1),
stage_id: i as i32,
mean_mw: 100.0,
std_mw: 0.0,
})
.collect();
let n_st = n_stages.max(1);
let mut bounds = ResolvedBounds::new(
1,
0,
0,
0,
0,
n_st,
HydroStageBounds {
min_storage_hm3: 0.0,
max_storage_hm3: 200.0,
min_turbined_m3s: 0.0,
max_turbined_m3s: 100.0,
min_outflow_m3s: 0.0,
max_outflow_m3s: None,
min_generation_mw: 0.0,
max_generation_mw: 250.0,
max_diversion_m3s: None,
filling_inflow_m3s: 0.0,
water_withdrawal_m3s,
},
ThermalStageBounds {
min_generation_mw: 0.0,
max_generation_mw: 0.0,
},
LineStageBounds {
direct_mw: 0.0,
reverse_mw: 0.0,
},
PumpingStageBounds {
min_flow_m3s: 0.0,
max_flow_m3s: 0.0,
},
ContractStageBounds {
min_mw: 0.0,
max_mw: 0.0,
price_per_mwh: 0.0,
},
);
for s in 0..n_st {
bounds.hydro_bounds_mut(0, s).water_withdrawal_m3s = water_withdrawal_m3s;
}
let mut penalties = ResolvedPenalties::new(
1,
1,
0,
0,
n_st,
HydroStagePenalties {
spillage_cost: 0.01,
diversion_cost: 0.0,
fpha_turbined_cost: 0.0,
storage_violation_below_cost: 0.0,
filling_target_violation_cost: 0.0,
turbined_violation_below_cost: 0.0,
outflow_violation_below_cost: 0.0,
outflow_violation_above_cost: 0.0,
generation_violation_below_cost: 0.0,
evaporation_violation_cost: 0.0,
water_withdrawal_violation_cost,
},
BusStagePenalties { excess_cost: 0.0 },
LineStagePenalties { exchange_cost: 0.0 },
NcsStagePenalties {
curtailment_cost: 0.0,
},
);
for s in 0..n_st {
penalties
.hydro_penalties_mut(0, s)
.water_withdrawal_violation_cost = water_withdrawal_violation_cost;
}
SystemBuilder::new()
.buses(vec![bus])
.hydros(vec![hydro])
.stages(stages)
.inflow_models(inflow_models)
.load_models(load_models)
.bounds(bounds)
.penalties(penalties)
.build()
.expect("one_hydro_system_with_withdrawal: valid")
}
#[test]
fn withdrawal_rhs_subtracted_from_water_balance() {
let withdrawal = 10.0_f64;
let system = one_hydro_system_with_withdrawal(1, 0, withdrawal, 0.0);
let result = build_stage_templates(
&system,
&no_penalty_config(),
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&default_production(&system),
&default_evaporation(&system),
)
.expect("build ok");
let t = &result.templates[0];
let row_water = 2_usize;
let total_hours = 744.0_f64;
let zeta = total_hours * 3_600.0 / 1_000_000.0;
let expected_rhs = zeta * (0.0 - withdrawal);
assert!(
(t.row_lower[row_water] - expected_rhs).abs() < 1e-12,
"water balance row_lower: expected {expected_rhs}, got {}",
t.row_lower[row_water]
);
assert!(
(t.row_upper[row_water] - expected_rhs).abs() < 1e-12,
"water balance row_upper: expected {expected_rhs}, got {}",
t.row_upper[row_water]
);
}
#[test]
fn withdrawal_zero_leaves_rhs_unchanged_from_base() {
let system_zero = one_hydro_system_with_withdrawal(1, 0, 0.0, 0.0);
let system_base = one_hydro_system(1, 0);
let result_zero = build_stage_templates(
&system_zero,
&no_penalty_config(),
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&default_production(&system_zero),
&default_evaporation(&system_zero),
)
.expect("zero-withdrawal build ok");
let result_base = build_stage_templates(
&system_base,
&no_penalty_config(),
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&default_production(&system_base),
&default_evaporation(&system_base),
)
.expect("base build ok");
let t_zero = &result_zero.templates[0];
let t_base = &result_base.templates[0];
let row_water = 1_usize;
assert!(
(t_zero.row_lower[row_water] - t_base.row_lower[row_water]).abs() < 1e-15,
"zero-withdrawal row_lower must equal base: {} vs {}",
t_zero.row_lower[row_water],
t_base.row_lower[row_water]
);
assert!(
(t_zero.row_upper[row_water] - t_base.row_upper[row_water]).abs() < 1e-15,
"zero-withdrawal row_upper must equal base: {} vs {}",
t_zero.row_upper[row_water],
t_base.row_upper[row_water]
);
}
#[test]
fn withdrawal_slack_matrix_entry_coefficient_is_minus_zeta() {
let system = one_hydro_system_with_withdrawal(1, 0, 5.0, 1000.0);
let result = build_stage_templates(
&system,
&no_penalty_config(),
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&default_production(&system),
&default_evaporation(&system),
)
.expect("build ok");
let t = &result.templates[0];
let col_w = t.num_cols - 1;
let row_water = 2_usize;
let total_hours = 744.0_f64;
let zeta = total_hours * 3_600.0 / 1_000_000.0;
let coeff = csc_entry(t, col_w, row_water).unwrap_or_else(|| {
panic!("withdrawal slack column {col_w} has no entry at water balance row {row_water}")
});
assert!(
(coeff - (-zeta)).abs() < 1e-12,
"withdrawal slack coefficient: expected {}, got {coeff}",
-zeta
);
let all_entries = entries_for_col(t, col_w);
assert_eq!(
all_entries.len(),
1,
"withdrawal slack column must have exactly 1 CSC entry, got {all_entries:?}"
);
}
#[test]
fn withdrawal_slack_objective_equals_cost_times_hours() {
let violation_cost = 1_000.0_f64;
let total_hours = 744.0_f64;
let system = one_hydro_system_with_withdrawal(1, 0, 5.0, violation_cost);
let result = build_stage_templates(
&system,
&no_penalty_config(),
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&default_production(&system),
&default_evaporation(&system),
)
.expect("build ok");
let t = &result.templates[0];
let col_w = t.num_cols - 1;
let expected_obj = violation_cost * total_hours / COST_SCALE_FACTOR;
assert!(
(t.objective[col_w] - expected_obj).abs() < 1e-12,
"withdrawal slack objective: expected {expected_obj}, got {}",
t.objective[col_w]
);
}
#[test]
fn withdrawal_slack_objective_zero_when_cost_is_zero() {
let system = one_hydro_system_with_withdrawal(1, 0, 0.0, 0.0);
let result = build_stage_templates(
&system,
&no_penalty_config(),
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&default_production(&system),
&default_evaporation(&system),
)
.expect("build ok");
let t = &result.templates[0];
let col_w = t.num_cols - 1;
assert!(
t.objective[col_w].abs() < 1e-15,
"withdrawal slack objective must be 0 when cost=0, got {}",
t.objective[col_w]
);
}
#[test]
fn withdrawal_slack_bounds_are_zero_to_infinity() {
let system = one_hydro_system_with_withdrawal(1, 0, 10.0, 5_000.0);
let result = build_stage_templates(
&system,
&no_penalty_config(),
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&default_production(&system),
&default_evaporation(&system),
)
.expect("build ok");
let t = &result.templates[0];
let col_w = t.num_cols - 1;
assert!(
(t.col_lower[col_w] - 0.0).abs() < 1e-15,
"withdrawal slack lower bound must be 0.0, got {}",
t.col_lower[col_w]
);
assert!(
t.col_upper[col_w].is_infinite() && t.col_upper[col_w] > 0.0,
"withdrawal slack upper bound must be +inf, got {}",
t.col_upper[col_w]
);
}
#[test]
#[allow(
clippy::cast_possible_truncation,
clippy::cast_possible_wrap,
clippy::too_many_lines
)]
fn two_hydro_withdrawal_slack_entries_per_hydro() {
use chrono::NaiveDate;
use cobre_core::entities::hydro::{Hydro, HydroGenerationModel, HydroPenalties};
use cobre_core::scenario::{InflowModel, LoadModel};
use cobre_core::temporal::{
Block, BlockMode, NoiseMethod, ScenarioSourceConfig, Stage, StageRiskConfig,
StageStateConfig,
};
let bus = Bus {
id: EntityId(1),
name: "B1".to_string(),
deficit_segments: vec![DeficitSegment {
depth_mw: None,
cost_per_mwh: 500.0,
}],
excess_cost: 0.0,
};
#[allow(clippy::cast_possible_wrap)]
let make_hydro = |id: i32| Hydro {
id: EntityId(id),
name: format!("H{id}"),
bus_id: EntityId(1),
downstream_id: None,
entry_stage_id: None,
exit_stage_id: None,
min_storage_hm3: 0.0,
max_storage_hm3: 200.0,
min_outflow_m3s: 0.0,
max_outflow_m3s: None,
generation_model: HydroGenerationModel::ConstantProductivity {
productivity_mw_per_m3s: 2.5,
},
min_turbined_m3s: 0.0,
max_turbined_m3s: 100.0,
min_generation_mw: 0.0,
max_generation_mw: 250.0,
tailrace: None,
hydraulic_losses: None,
efficiency: None,
evaporation_coefficients_mm: None,
evaporation_reference_volumes_hm3: None,
diversion: None,
filling: None,
penalties: HydroPenalties {
spillage_cost: 0.01,
diversion_cost: 0.0,
fpha_turbined_cost: 0.0,
storage_violation_below_cost: 0.0,
filling_target_violation_cost: 0.0,
turbined_violation_below_cost: 0.0,
outflow_violation_below_cost: 0.0,
outflow_violation_above_cost: 0.0,
generation_violation_below_cost: 0.0,
evaporation_violation_cost: 0.0,
water_withdrawal_violation_cost: 2_000.0,
},
};
let stages = vec![Stage {
index: 0,
id: 0,
start_date: NaiveDate::from_ymd_opt(2024, 1, 1).unwrap(),
end_date: NaiveDate::from_ymd_opt(2024, 2, 1).unwrap(),
season_id: None,
blocks: vec![Block {
index: 0,
name: "S".to_string(),
duration_hours: 744.0,
}],
block_mode: BlockMode::Parallel,
state_config: StageStateConfig {
storage: true,
inflow_lags: false,
},
risk_config: StageRiskConfig::Expectation,
scenario_config: ScenarioSourceConfig {
branching_factor: 1,
noise_method: NoiseMethod::Saa,
},
}];
let inflow_models = vec![
InflowModel {
hydro_id: EntityId(2),
stage_id: 0,
mean_m3s: 80.0,
std_m3s: 20.0,
ar_coefficients: vec![],
residual_std_ratio: 1.0,
},
InflowModel {
hydro_id: EntityId(3),
stage_id: 0,
mean_m3s: 50.0,
std_m3s: 10.0,
ar_coefficients: vec![],
residual_std_ratio: 1.0,
},
];
let load_models = vec![LoadModel {
bus_id: EntityId(1),
stage_id: 0,
mean_mw: 100.0,
std_mw: 0.0,
}];
let hydro_bounds_default = HydroStageBounds {
min_storage_hm3: 0.0,
max_storage_hm3: 200.0,
min_turbined_m3s: 0.0,
max_turbined_m3s: 100.0,
min_outflow_m3s: 0.0,
max_outflow_m3s: None,
min_generation_mw: 0.0,
max_generation_mw: 250.0,
max_diversion_m3s: None,
filling_inflow_m3s: 0.0,
water_withdrawal_m3s: 0.0,
};
let mut bounds = ResolvedBounds::new(
2,
0,
0,
0,
0,
1,
hydro_bounds_default,
ThermalStageBounds {
min_generation_mw: 0.0,
max_generation_mw: 0.0,
},
LineStageBounds {
direct_mw: 0.0,
reverse_mw: 0.0,
},
PumpingStageBounds {
min_flow_m3s: 0.0,
max_flow_m3s: 0.0,
},
ContractStageBounds {
min_mw: 0.0,
max_mw: 0.0,
price_per_mwh: 0.0,
},
);
bounds.hydro_bounds_mut(0, 0).water_withdrawal_m3s = 8.0;
bounds.hydro_bounds_mut(1, 0).water_withdrawal_m3s = 12.0;
let hydro_penalties_default = HydroStagePenalties {
spillage_cost: 0.01,
diversion_cost: 0.0,
fpha_turbined_cost: 0.0,
storage_violation_below_cost: 0.0,
filling_target_violation_cost: 0.0,
turbined_violation_below_cost: 0.0,
outflow_violation_below_cost: 0.0,
outflow_violation_above_cost: 0.0,
generation_violation_below_cost: 0.0,
evaporation_violation_cost: 0.0,
water_withdrawal_violation_cost: 2_000.0,
};
let penalties = ResolvedPenalties::new(
2,
1,
0,
0,
1,
hydro_penalties_default,
BusStagePenalties { excess_cost: 0.0 },
LineStagePenalties { exchange_cost: 0.0 },
NcsStagePenalties {
curtailment_cost: 0.0,
},
);
let system = SystemBuilder::new()
.buses(vec![bus])
.hydros(vec![make_hydro(2), make_hydro(3)])
.stages(stages)
.inflow_models(inflow_models)
.load_models(load_models)
.bounds(bounds)
.penalties(penalties)
.build()
.expect("two_hydro_with_withdrawal: valid");
let result = build_stage_templates(
&system,
&no_penalty_config(),
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&default_production(&system),
&default_evaporation(&system),
)
.expect("build ok");
let t = &result.templates[0];
let col_w0 = t.num_cols - 2;
let col_w1 = t.num_cols - 1;
let row_w0 = 4_usize; let row_w1 = 5_usize;
let total_hours = 744.0_f64;
let zeta = total_hours * 3_600.0 / 1_000_000.0;
let coeff_w0 = csc_entry(t, col_w0, row_w0).unwrap_or_else(|| {
panic!("withdrawal slack col {col_w0} has no entry at water balance row {row_w0}")
});
assert!(
(coeff_w0 - (-zeta)).abs() < 1e-12,
"hydro-0 withdrawal slack coeff: expected {}, got {coeff_w0}",
-zeta
);
let coeff_w1 = csc_entry(t, col_w1, row_w1).unwrap_or_else(|| {
panic!("withdrawal slack col {col_w1} has no entry at water balance row {row_w1}")
});
assert!(
(coeff_w1 - (-zeta)).abs() < 1e-12,
"hydro-1 withdrawal slack coeff: expected {}, got {coeff_w1}",
-zeta
);
assert!(
csc_entry(t, col_w0, row_w1).is_none(),
"hydro-0 withdrawal slack must not appear in hydro-1 water balance row"
);
assert!(
csc_entry(t, col_w1, row_w0).is_none(),
"hydro-1 withdrawal slack must not appear in hydro-0 water balance row"
);
}
#[test]
#[allow(clippy::too_many_lines)]
fn three_hydro_num_cols_includes_three_withdrawal_slacks() {
use chrono::NaiveDate;
use cobre_core::entities::hydro::{Hydro, HydroGenerationModel, HydroPenalties};
use cobre_core::scenario::{InflowModel, LoadModel};
use cobre_core::temporal::{
Block, BlockMode, NoiseMethod, ScenarioSourceConfig, Stage, StageRiskConfig,
StageStateConfig,
};
let bus = Bus {
id: EntityId(1),
name: "B1".to_string(),
deficit_segments: vec![DeficitSegment {
depth_mw: None,
cost_per_mwh: 500.0,
}],
excess_cost: 0.0,
};
#[allow(clippy::cast_possible_wrap)]
let make_hydro = |id: i32| Hydro {
id: EntityId(id),
name: format!("H{id}"),
bus_id: EntityId(1),
downstream_id: None,
entry_stage_id: None,
exit_stage_id: None,
min_storage_hm3: 0.0,
max_storage_hm3: 200.0,
min_outflow_m3s: 0.0,
max_outflow_m3s: None,
generation_model: HydroGenerationModel::ConstantProductivity {
productivity_mw_per_m3s: 2.5,
},
min_turbined_m3s: 0.0,
max_turbined_m3s: 100.0,
min_generation_mw: 0.0,
max_generation_mw: 250.0,
tailrace: None,
hydraulic_losses: None,
efficiency: None,
evaporation_coefficients_mm: None,
evaporation_reference_volumes_hm3: None,
diversion: None,
filling: None,
penalties: HydroPenalties {
spillage_cost: 0.01,
diversion_cost: 0.0,
fpha_turbined_cost: 0.0,
storage_violation_below_cost: 0.0,
filling_target_violation_cost: 0.0,
turbined_violation_below_cost: 0.0,
outflow_violation_below_cost: 0.0,
outflow_violation_above_cost: 0.0,
generation_violation_below_cost: 0.0,
evaporation_violation_cost: 0.0,
water_withdrawal_violation_cost: 1_000.0,
},
};
let stages = vec![Stage {
index: 0,
id: 0,
start_date: NaiveDate::from_ymd_opt(2024, 1, 1).unwrap(),
end_date: NaiveDate::from_ymd_opt(2024, 2, 1).unwrap(),
season_id: None,
blocks: vec![Block {
index: 0,
name: "S".to_string(),
duration_hours: 744.0,
}],
block_mode: BlockMode::Parallel,
state_config: StageStateConfig {
storage: true,
inflow_lags: false,
},
risk_config: StageRiskConfig::Expectation,
scenario_config: ScenarioSourceConfig {
branching_factor: 1,
noise_method: NoiseMethod::Saa,
},
}];
let inflow_models: Vec<InflowModel> = [1, 2, 3]
.iter()
.map(|&hid| InflowModel {
hydro_id: EntityId(hid),
stage_id: 0,
mean_m3s: 50.0,
std_m3s: 10.0,
ar_coefficients: vec![],
residual_std_ratio: 1.0,
})
.collect();
let load_models = vec![LoadModel {
bus_id: EntityId(1),
stage_id: 0,
mean_mw: 100.0,
std_mw: 0.0,
}];
let bounds = ResolvedBounds::new(
3,
0,
0,
0,
0,
1,
HydroStageBounds {
min_storage_hm3: 0.0,
max_storage_hm3: 200.0,
min_turbined_m3s: 0.0,
max_turbined_m3s: 100.0,
min_outflow_m3s: 0.0,
max_outflow_m3s: None,
min_generation_mw: 0.0,
max_generation_mw: 250.0,
max_diversion_m3s: None,
filling_inflow_m3s: 0.0,
water_withdrawal_m3s: 5.0,
},
ThermalStageBounds {
min_generation_mw: 0.0,
max_generation_mw: 0.0,
},
LineStageBounds {
direct_mw: 0.0,
reverse_mw: 0.0,
},
PumpingStageBounds {
min_flow_m3s: 0.0,
max_flow_m3s: 0.0,
},
ContractStageBounds {
min_mw: 0.0,
max_mw: 0.0,
price_per_mwh: 0.0,
},
);
let penalties = ResolvedPenalties::new(
3,
1,
0,
0,
1,
HydroStagePenalties {
spillage_cost: 0.01,
diversion_cost: 0.0,
fpha_turbined_cost: 0.0,
storage_violation_below_cost: 0.0,
filling_target_violation_cost: 0.0,
turbined_violation_below_cost: 0.0,
outflow_violation_below_cost: 0.0,
outflow_violation_above_cost: 0.0,
generation_violation_below_cost: 0.0,
evaporation_violation_cost: 0.0,
water_withdrawal_violation_cost: 1_000.0,
},
BusStagePenalties { excess_cost: 0.0 },
LineStagePenalties { exchange_cost: 0.0 },
NcsStagePenalties {
curtailment_cost: 0.0,
},
);
let system = SystemBuilder::new()
.buses(vec![bus])
.hydros(vec![make_hydro(1), make_hydro(2), make_hydro(3)])
.stages(stages)
.inflow_models(inflow_models)
.load_models(load_models)
.bounds(bounds)
.penalties(penalties)
.build()
.expect("three_hydro_system: valid");
let result = build_stage_templates(
&system,
&no_penalty_config(),
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&default_production(&system),
&default_evaporation(&system),
)
.expect("build ok");
let t = &result.templates[0];
let expected_cols = 24_usize;
assert_eq!(
t.num_cols, expected_cols,
"3-hydro system: num_cols should be {expected_cols}, got {}",
t.num_cols
);
let n_h = t.n_hydro;
assert_eq!(
t.col_upper[t.num_cols - n_h],
f64::INFINITY,
"withdrawal slack column for hydro 0 should be unbounded above"
);
assert_eq!(
t.col_upper[t.num_cols - n_h + 1],
f64::INFINITY,
"withdrawal slack column for hydro 1 should be unbounded above"
);
assert_eq!(
t.col_upper[t.num_cols - n_h + 2],
f64::INFINITY,
"withdrawal slack column for hydro 2 should be unbounded above"
);
}
#[allow(clippy::cast_possible_wrap)]
fn one_bus_system_n_blks(n_blks: usize) -> cobre_core::System {
use chrono::NaiveDate;
use cobre_core::scenario::LoadModel;
use cobre_core::temporal::{
Block, BlockMode, NoiseMethod, ScenarioSourceConfig, StageRiskConfig, StageStateConfig,
};
let bus = Bus {
id: EntityId(1),
name: "B1".to_string(),
deficit_segments: vec![DeficitSegment {
depth_mw: None,
cost_per_mwh: 500.0,
}],
excess_cost: 0.0,
};
let blocks: Vec<Block> = (0..n_blks)
.map(|i| Block {
index: i,
name: format!("BLK{i}"),
duration_hours: 720.0,
})
.collect();
let stage = cobre_core::temporal::Stage {
index: 0,
id: 0_i32,
start_date: NaiveDate::from_ymd_opt(2024, 1, 1).unwrap(),
end_date: NaiveDate::from_ymd_opt(2024, 2, 1).unwrap(),
season_id: None,
blocks,
block_mode: BlockMode::Parallel,
state_config: StageStateConfig {
storage: false,
inflow_lags: false,
},
risk_config: StageRiskConfig::Expectation,
scenario_config: ScenarioSourceConfig {
branching_factor: 1,
noise_method: NoiseMethod::Saa,
},
};
let load_models: Vec<LoadModel> = vec![LoadModel {
bus_id: EntityId(1),
stage_id: 0,
mean_mw: 100.0,
std_mw: 0.0,
}];
let bounds = ResolvedBounds::new(
0,
0,
0,
0,
0,
1,
default_hydro_bounds(),
ThermalStageBounds {
min_generation_mw: 0.0,
max_generation_mw: 0.0,
},
LineStageBounds {
direct_mw: 0.0,
reverse_mw: 0.0,
},
PumpingStageBounds {
min_flow_m3s: 0.0,
max_flow_m3s: 0.0,
},
ContractStageBounds {
min_mw: 0.0,
max_mw: 0.0,
price_per_mwh: 0.0,
},
);
let penalties = ResolvedPenalties::new(
0,
1,
0,
0,
1,
default_hydro_penalties(),
BusStagePenalties { excess_cost: 0.0 },
LineStagePenalties { exchange_cost: 0.0 },
NcsStagePenalties {
curtailment_cost: 0.0,
},
);
SystemBuilder::new()
.buses(vec![bus])
.stages(vec![stage])
.load_models(load_models)
.bounds(bounds)
.penalties(penalties)
.build()
.expect("one_bus_system_n_blks: valid")
}
fn make_constraint(
id: i32,
sense: cobre_core::ConstraintSense,
slack_enabled: bool,
) -> cobre_core::GenericConstraint {
use cobre_core::{ConstraintExpression, GenericConstraint, SlackConfig};
GenericConstraint {
id: EntityId(id),
name: format!("gc_{id}"),
description: None,
expression: ConstraintExpression { terms: vec![] },
sense,
slack: SlackConfig {
enabled: slack_enabled,
penalty: if slack_enabled { Some(5000.0) } else { None },
},
}
}
fn build_templates_for(system: &cobre_core::System) -> Vec<cobre_solver::StageTemplate> {
let production = default_production(system);
let evaporation = default_evaporation(system);
build_stage_templates(
system,
&no_penalty_config(),
&PrecomputedPar::default(),
&PrecomputedNormal::default(),
&production,
&evaporation,
)
.expect("build_templates_for: valid")
.templates
}
#[test]
fn generic_constraints_zero_does_not_change_layout() {
let system = one_bus_system_n_blks(1);
let templates = build_templates_for(&system);
let t = &templates[0];
assert!(t.num_cols > 0);
assert!(t.num_rows > 0);
let templates2 = build_templates_for(&system);
assert_eq!(
templates2[0].num_cols, t.num_cols,
"second build must not change num_cols"
);
assert_eq!(
templates2[0].num_rows, t.num_rows,
"second build must not change num_rows"
);
}
#[test]
fn generic_constraint_no_slack_block_id_none_3_blocks() {
use cobre_core::{ConstraintSense, ResolvedGenericConstraintBounds};
use std::collections::HashMap;
let n_blks = 3_usize;
let baseline_system = one_bus_system_n_blks(n_blks);
let baseline_rows = build_templates_for(&baseline_system)[0].num_rows;
let baseline_cols = build_templates_for(&baseline_system)[0].num_cols;
let id_map: HashMap<i32, usize> = [(10_i32, 0)].into_iter().collect();
let rows = vec![(10_i32, 0_i32, None::<i32>, 500.0_f64)];
let generic_bounds = ResolvedGenericConstraintBounds::new(&id_map, rows.into_iter());
let constraint = make_constraint(10, ConstraintSense::LessEqual, false);
let system = one_bus_system_n_blks_with_generic(n_blks, vec![constraint], generic_bounds);
let t_rows = build_templates_for(&system)[0].num_rows;
let t_cols = build_templates_for(&system)[0].num_cols;
assert_eq!(
t_rows,
baseline_rows + n_blks,
"num_rows must increase by n_blks={n_blks} (one row per block, no slack)"
);
assert_eq!(
t_cols, baseline_cols,
"num_cols must be unchanged (no slack columns)"
);
}
#[test]
fn generic_constraint_le_slack_enabled_2_blocks() {
use cobre_core::{ConstraintSense, ResolvedGenericConstraintBounds};
use std::collections::HashMap;
let n_blks = 2_usize;
let baseline_system = one_bus_system_n_blks(n_blks);
let baseline_rows = build_templates_for(&baseline_system)[0].num_rows;
let baseline_cols = build_templates_for(&baseline_system)[0].num_cols;
let id_map: HashMap<i32, usize> = [(20_i32, 0)].into_iter().collect();
let rows = vec![(20_i32, 0_i32, None::<i32>, 300.0_f64)];
let generic_bounds = ResolvedGenericConstraintBounds::new(&id_map, rows.into_iter());
let constraint = make_constraint(20, ConstraintSense::LessEqual, true);
let system = one_bus_system_n_blks_with_generic(n_blks, vec![constraint], generic_bounds);
let t_rows = build_templates_for(&system)[0].num_rows;
let t_cols = build_templates_for(&system)[0].num_cols;
assert_eq!(
t_rows,
baseline_rows + 2,
"num_rows must increase by 2 (one row per block)"
);
assert_eq!(
t_cols,
baseline_cols + 2,
"num_cols must increase by 2 (one slack per row for `<=`)"
);
}
#[test]
fn generic_constraint_equal_sense_two_slacks_per_row() {
use cobre_core::{ConstraintSense, ResolvedGenericConstraintBounds};
use std::collections::HashMap;
let n_blks = 2_usize;
let baseline_system = one_bus_system_n_blks(n_blks);
let baseline_rows = build_templates_for(&baseline_system)[0].num_rows;
let baseline_cols = build_templates_for(&baseline_system)[0].num_cols;
let id_map: HashMap<i32, usize> = [(30_i32, 0)].into_iter().collect();
let rows = vec![(30_i32, 0_i32, None::<i32>, 100.0_f64)];
let generic_bounds = ResolvedGenericConstraintBounds::new(&id_map, rows.into_iter());
let constraint = make_constraint(30, ConstraintSense::Equal, true);
let system = one_bus_system_n_blks_with_generic(n_blks, vec![constraint], generic_bounds);
let t_rows = build_templates_for(&system)[0].num_rows;
let t_cols = build_templates_for(&system)[0].num_cols;
assert_eq!(
t_rows,
baseline_rows + 2,
"num_rows must increase by 2 (one row per block)"
);
assert_eq!(
t_cols,
baseline_cols + 4,
"num_cols must increase by 4 (two slacks per row for `==`)"
);
}
#[test]
fn generic_constraint_specific_block_id_generates_one_row() {
use cobre_core::{ConstraintSense, ResolvedGenericConstraintBounds};
use std::collections::HashMap;
let n_blks = 3_usize;
let baseline_system = one_bus_system_n_blks(n_blks);
let baseline_rows = build_templates_for(&baseline_system)[0].num_rows;
let baseline_cols = build_templates_for(&baseline_system)[0].num_cols;
let id_map: HashMap<i32, usize> = [(40_i32, 0)].into_iter().collect();
let rows = vec![(40_i32, 0_i32, Some(1_i32), 200.0_f64)];
let generic_bounds = ResolvedGenericConstraintBounds::new(&id_map, rows.into_iter());
let constraint = make_constraint(40, ConstraintSense::LessEqual, false);
let system = one_bus_system_n_blks_with_generic(n_blks, vec![constraint], generic_bounds);
let t_rows = build_templates_for(&system)[0].num_rows;
let t_cols = build_templates_for(&system)[0].num_cols;
assert_eq!(
t_rows,
baseline_rows + 1,
"num_rows must increase by exactly 1 (only the specified block)"
);
assert_eq!(
t_cols, baseline_cols,
"num_cols must be unchanged (no slack columns)"
);
}
#[test]
fn generic_constraint_inactive_does_not_contribute_rows() {
use cobre_core::{ConstraintSense, ResolvedGenericConstraintBounds};
use std::collections::HashMap;
let n_blks = 2_usize;
let baseline_system = one_bus_system_n_blks(n_blks);
let baseline_rows = build_templates_for(&baseline_system)[0].num_rows;
let id_map: HashMap<i32, usize> = [(50_i32, 0), (51_i32, 1)].into_iter().collect();
let rows = vec![(50_i32, 0_i32, None::<i32>, 400.0_f64)];
let generic_bounds = ResolvedGenericConstraintBounds::new(&id_map, rows.into_iter());
let c_active = make_constraint(50, ConstraintSense::LessEqual, false);
let c_inactive = make_constraint(51, ConstraintSense::LessEqual, false);
let system =
one_bus_system_n_blks_with_generic(n_blks, vec![c_active, c_inactive], generic_bounds);
let t_rows = build_templates_for(&system)[0].num_rows;
assert_eq!(
t_rows,
baseline_rows + n_blks,
"only the active constraint must contribute rows"
);
}
#[test]
fn stage_indexer_generic_fields_empty_from_new() {
let idx = crate::indexer::StageIndexer::new(3, 2);
assert!(
idx.generic_constraint_rows.is_empty(),
"generic_constraint_rows must be empty from new()"
);
assert!(
idx.generic_constraint_slack.is_empty(),
"generic_constraint_slack must be empty from new()"
);
assert_eq!(
idx.n_generic_constraints_active, 0,
"n_generic_constraints_active must be 0 from new()"
);
}
#[allow(clippy::cast_possible_wrap)]
fn one_bus_system_n_blks_with_generic(
n_blks: usize,
constraints: Vec<cobre_core::GenericConstraint>,
bounds: cobre_core::ResolvedGenericConstraintBounds,
) -> cobre_core::System {
use chrono::NaiveDate;
use cobre_core::scenario::LoadModel;
use cobre_core::temporal::{
Block, BlockMode, NoiseMethod, ScenarioSourceConfig, Stage, StageRiskConfig,
StageStateConfig,
};
let bus = Bus {
id: EntityId(1),
name: "B1".to_string(),
deficit_segments: vec![DeficitSegment {
depth_mw: None,
cost_per_mwh: 500.0,
}],
excess_cost: 0.0,
};
let blks: Vec<Block> = (0..n_blks)
.map(|i| Block {
index: i,
name: format!("BLK{i}"),
duration_hours: 720.0,
})
.collect();
let stage = Stage {
index: 0,
id: 0_i32,
start_date: NaiveDate::from_ymd_opt(2024, 1, 1).unwrap(),
end_date: NaiveDate::from_ymd_opt(2024, 2, 1).unwrap(),
season_id: None,
blocks: blks,
block_mode: BlockMode::Parallel,
state_config: StageStateConfig {
storage: false,
inflow_lags: false,
},
risk_config: StageRiskConfig::Expectation,
scenario_config: ScenarioSourceConfig {
branching_factor: 1,
noise_method: NoiseMethod::Saa,
},
};
let load_models: Vec<LoadModel> = vec![LoadModel {
bus_id: EntityId(1),
stage_id: 0,
mean_mw: 100.0,
std_mw: 0.0,
}];
let resolved_bounds = ResolvedBounds::new(
0,
0,
0,
0,
0,
1,
default_hydro_bounds(),
ThermalStageBounds {
min_generation_mw: 0.0,
max_generation_mw: 0.0,
},
LineStageBounds {
direct_mw: 0.0,
reverse_mw: 0.0,
},
PumpingStageBounds {
min_flow_m3s: 0.0,
max_flow_m3s: 0.0,
},
ContractStageBounds {
min_mw: 0.0,
max_mw: 0.0,
price_per_mwh: 0.0,
},
);
let penalties = ResolvedPenalties::new(
0,
1,
0,
0,
1,
default_hydro_penalties(),
BusStagePenalties { excess_cost: 0.0 },
LineStagePenalties { exchange_cost: 0.0 },
NcsStagePenalties {
curtailment_cost: 0.0,
},
);
SystemBuilder::new()
.buses(vec![bus])
.stages(vec![stage])
.load_models(load_models)
.bounds(resolved_bounds)
.penalties(penalties)
.generic_constraints(constraints)
.resolved_generic_bounds(bounds)
.build()
.expect("one_bus_system_n_blks_with_generic: valid")
}
fn csc_entries_at(t: &cobre_solver::StageTemplate, col: usize, row: usize) -> Vec<f64> {
let start = t.col_starts[col] as usize;
let end = t.col_starts[col + 1] as usize;
t.row_indices[start..end]
.iter()
.zip(t.values[start..end].iter())
.filter_map(|(&r, &v)| if r as usize == row { Some(v) } else { None })
.collect()
}
#[test]
fn generic_constraint_thermal_le_row_bounds_and_csc_entry() {
use cobre_core::ResolvedGenericConstraintBounds;
use cobre_core::{
ConstraintExpression, ConstraintSense, GenericConstraint, LinearTerm, SlackConfig,
VariableRef,
};
use std::collections::HashMap;
let thermal_entity_id = EntityId(2);
let constraint = GenericConstraint {
id: EntityId(10),
name: "gc_thermal_le".to_string(),
description: None,
expression: ConstraintExpression {
terms: vec![LinearTerm {
coefficient: 1.0,
variable: VariableRef::ThermalGeneration {
thermal_id: thermal_entity_id,
block_id: None,
},
}],
},
sense: ConstraintSense::LessEqual,
slack: SlackConfig {
enabled: false,
penalty: None,
},
};
let id_map: HashMap<i32, usize> = [(10_i32, 0)].into_iter().collect();
let rows = vec![(10_i32, 0_i32, None::<i32>, 50.0_f64)];
let generic_bounds = ResolvedGenericConstraintBounds::new(&id_map, rows.into_iter());
let system =
one_bus_one_thermal_system(thermal_entity_id, vec![constraint], generic_bounds);
let t = &build_templates_for(&system)[0];
let thermal_col = 1_usize;
let generic_row = 1_usize;
assert!(
t.row_lower[generic_row].is_infinite() && t.row_lower[generic_row] < 0.0,
"row_lower must be -INF for <= constraint, got {}",
t.row_lower[generic_row]
);
assert!(
(t.row_upper[generic_row] - 50.0).abs() < f64::EPSILON,
"row_upper must be 50.0, got {}",
t.row_upper[generic_row]
);
let entries = csc_entries_at(t, thermal_col, generic_row);
assert!(
!entries.is_empty(),
"no CSC entry found at (col={thermal_col}, row={generic_row})"
);
let total: f64 = entries.iter().sum();
assert!(
(total - 1.0).abs() < f64::EPSILON,
"expected +1.0 total at thermal col / generic row, got {total}"
);
}
#[test]
fn generic_constraint_thermal_le_slack_column_and_csc_entry() {
use cobre_core::ResolvedGenericConstraintBounds;
use cobre_core::{
ConstraintExpression, ConstraintSense, GenericConstraint, LinearTerm, SlackConfig,
VariableRef,
};
use std::collections::HashMap;
let thermal_entity_id = EntityId(2);
let block_hours = 744.0_f64;
let penalty = 5000.0_f64;
let constraint = GenericConstraint {
id: EntityId(20),
name: "gc_thermal_le_slack".to_string(),
description: None,
expression: ConstraintExpression {
terms: vec![LinearTerm {
coefficient: 1.0,
variable: VariableRef::ThermalGeneration {
thermal_id: thermal_entity_id,
block_id: None,
},
}],
},
sense: ConstraintSense::LessEqual,
slack: SlackConfig {
enabled: true,
penalty: Some(penalty),
},
};
let id_map: HashMap<i32, usize> = [(20_i32, 0)].into_iter().collect();
let rows = vec![(20_i32, 0_i32, None::<i32>, 50.0_f64)];
let generic_bounds = ResolvedGenericConstraintBounds::new(&id_map, rows.into_iter());
let system =
one_bus_one_thermal_system(thermal_entity_id, vec![constraint], generic_bounds);
let t = &build_templates_for(&system)[0];
let slack_col = 4_usize;
let generic_row = 1_usize;
assert!(
t.col_lower[slack_col].abs() < f64::EPSILON,
"slack col_lower must be 0.0, got {}",
t.col_lower[slack_col]
);
assert!(
t.col_upper[slack_col].is_infinite() && t.col_upper[slack_col] > 0.0,
"slack col_upper must be +INF, got {}",
t.col_upper[slack_col]
);
let expected_obj = penalty * block_hours / COST_SCALE_FACTOR;
assert!(
(t.objective[slack_col] - expected_obj).abs() < 1e-12,
"slack objective must be {expected_obj}, got {}",
t.objective[slack_col]
);
let entries = csc_entries_at(t, slack_col, generic_row);
assert!(
!entries.is_empty(),
"no CSC entry found at (col={slack_col}, row={generic_row})"
);
let total: f64 = entries.iter().sum();
assert!(
(total - (-1.0)).abs() < f64::EPSILON,
"expected -1.0 at slack col / generic row for <= sense, got {total}"
);
}
#[test]
fn generic_constraint_thermal_ge_row_bounds() {
use cobre_core::ResolvedGenericConstraintBounds;
use cobre_core::{
ConstraintExpression, ConstraintSense, GenericConstraint, LinearTerm, SlackConfig,
VariableRef,
};
use std::collections::HashMap;
let thermal_entity_id = EntityId(2);
let constraint = GenericConstraint {
id: EntityId(30),
name: "gc_thermal_ge".to_string(),
description: None,
expression: ConstraintExpression {
terms: vec![LinearTerm {
coefficient: 1.0,
variable: VariableRef::ThermalGeneration {
thermal_id: thermal_entity_id,
block_id: None,
},
}],
},
sense: ConstraintSense::GreaterEqual,
slack: SlackConfig {
enabled: false,
penalty: None,
},
};
let id_map: HashMap<i32, usize> = [(30_i32, 0)].into_iter().collect();
let rows = vec![(30_i32, 0_i32, None::<i32>, 10.0_f64)];
let generic_bounds = ResolvedGenericConstraintBounds::new(&id_map, rows.into_iter());
let system =
one_bus_one_thermal_system(thermal_entity_id, vec![constraint], generic_bounds);
let t = &build_templates_for(&system)[0];
let generic_row = 1_usize;
assert!(
(t.row_lower[generic_row] - 10.0).abs() < f64::EPSILON,
"row_lower must be 10.0 for >= constraint, got {}",
t.row_lower[generic_row]
);
assert!(
t.row_upper[generic_row].is_infinite() && t.row_upper[generic_row] > 0.0,
"row_upper must be +INF for >= constraint, got {}",
t.row_upper[generic_row]
);
}
#[test]
fn generic_constraint_thermal_equal_two_slacks() {
use cobre_core::ResolvedGenericConstraintBounds;
use cobre_core::{ConstraintExpression, ConstraintSense, GenericConstraint, SlackConfig};
use std::collections::HashMap;
let thermal_entity_id = EntityId(2);
let penalty = 5000.0_f64;
let block_hours = 744.0_f64;
let constraint = GenericConstraint {
id: EntityId(40),
name: "gc_thermal_eq_slack".to_string(),
description: None,
expression: ConstraintExpression { terms: vec![] },
sense: ConstraintSense::Equal,
slack: SlackConfig {
enabled: true,
penalty: Some(penalty),
},
};
let id_map: HashMap<i32, usize> = [(40_i32, 0)].into_iter().collect();
let rows = vec![(40_i32, 0_i32, None::<i32>, 80.0_f64)];
let generic_bounds = ResolvedGenericConstraintBounds::new(&id_map, rows.into_iter());
let system =
one_bus_one_thermal_system(thermal_entity_id, vec![constraint], generic_bounds);
let t = &build_templates_for(&system)[0];
let slack_plus_col = 4_usize;
let slack_minus_col = 5_usize;
let generic_row = 1_usize;
assert!(
(t.row_lower[generic_row] - 80.0).abs() < f64::EPSILON,
"row_lower must be 80.0 for == constraint, got {}",
t.row_lower[generic_row]
);
assert!(
(t.row_upper[generic_row] - 80.0).abs() < f64::EPSILON,
"row_upper must be 80.0 for == constraint, got {}",
t.row_upper[generic_row]
);
assert_eq!(t.num_cols, 6, "num_cols must be 6 with 2 slack columns");
assert!(
t.col_upper[slack_plus_col].is_infinite() && t.col_upper[slack_plus_col] > 0.0,
"plus slack col_upper must be +INF"
);
let expected_obj = penalty * block_hours / COST_SCALE_FACTOR;
assert!(
(t.objective[slack_plus_col] - expected_obj).abs() < 1e-12,
"plus slack objective must be {expected_obj}, got {}",
t.objective[slack_plus_col]
);
let plus_entries = csc_entries_at(t, slack_plus_col, generic_row);
assert!(
!plus_entries.is_empty(),
"no CSC entry at plus slack col / generic row"
);
let plus_total: f64 = plus_entries.iter().sum();
assert!(
(plus_total - 1.0).abs() < f64::EPSILON,
"plus slack CSC must be +1.0 for == sense, got {plus_total}"
);
assert!(
t.col_upper[slack_minus_col].is_infinite() && t.col_upper[slack_minus_col] > 0.0,
"minus slack col_upper must be +INF"
);
assert!(
(t.objective[slack_minus_col] - expected_obj).abs() < 1e-12,
"minus slack objective must be {expected_obj}"
);
let minus_entries = csc_entries_at(t, slack_minus_col, generic_row);
assert!(
!minus_entries.is_empty(),
"no CSC entry at minus slack col / generic row"
);
let minus_total: f64 = minus_entries.iter().sum();
assert!(
(minus_total - (-1.0)).abs() < f64::EPSILON,
"minus slack CSC must be -1.0 for == sense, got {minus_total}"
);
}
#[test]
#[allow(clippy::cast_possible_wrap)]
fn generic_constraint_two_hydros_sum_csc_entries() {
use chrono::NaiveDate;
use cobre_core::ResolvedGenericConstraintBounds;
use cobre_core::entities::hydro::{Hydro, HydroGenerationModel, HydroPenalties};
use cobre_core::scenario::{InflowModel, LoadModel};
use cobre_core::temporal::{
Block, BlockMode, NoiseMethod, ScenarioSourceConfig, Stage, StageRiskConfig,
StageStateConfig,
};
use cobre_core::{
ConstraintExpression, ConstraintSense, GenericConstraint, LinearTerm, SlackConfig,
VariableRef,
};
use std::collections::HashMap;
let h1_id = EntityId(5);
let h2_id = EntityId(10);
let prod_h1 = 2.5_f64;
let prod_h2 = 3.0_f64;
let make_hydro = |id: EntityId, prod: f64| Hydro {
id,
name: format!("H{}", id.0),
bus_id: EntityId(1),
downstream_id: None,
entry_stage_id: None,
exit_stage_id: None,
min_storage_hm3: 0.0,
max_storage_hm3: 200.0,
min_outflow_m3s: 0.0,
max_outflow_m3s: None,
generation_model: HydroGenerationModel::ConstantProductivity {
productivity_mw_per_m3s: prod,
},
min_turbined_m3s: 0.0,
max_turbined_m3s: 100.0,
min_generation_mw: 0.0,
max_generation_mw: 250.0,
tailrace: None,
hydraulic_losses: None,
efficiency: None,
evaporation_coefficients_mm: None,
evaporation_reference_volumes_hm3: None,
diversion: None,
filling: None,
penalties: HydroPenalties {
spillage_cost: 0.01,
diversion_cost: 0.0,
fpha_turbined_cost: 0.0,
storage_violation_below_cost: 0.0,
filling_target_violation_cost: 0.0,
turbined_violation_below_cost: 0.0,
outflow_violation_below_cost: 0.0,
outflow_violation_above_cost: 0.0,
generation_violation_below_cost: 0.0,
evaporation_violation_cost: 0.0,
water_withdrawal_violation_cost: 0.0,
},
};
let hydros = vec![make_hydro(h1_id, prod_h1), make_hydro(h2_id, prod_h2)];
let bus = Bus {
id: EntityId(1),
name: "B1".to_string(),
deficit_segments: vec![DeficitSegment {
depth_mw: None,
cost_per_mwh: 500.0,
}],
excess_cost: 0.0,
};
let stage = Stage {
index: 0,
id: 0_i32,
start_date: NaiveDate::from_ymd_opt(2024, 1, 1).unwrap(),
end_date: NaiveDate::from_ymd_opt(2024, 2, 1).unwrap(),
season_id: None,
blocks: vec![Block {
index: 0,
name: "BLK0".to_string(),
duration_hours: 744.0,
}],
block_mode: BlockMode::Parallel,
state_config: StageStateConfig {
storage: true,
inflow_lags: false,
},
risk_config: StageRiskConfig::Expectation,
scenario_config: ScenarioSourceConfig {
branching_factor: 1,
noise_method: NoiseMethod::Saa,
},
};
let inflow_models = vec![
InflowModel {
hydro_id: h1_id,
stage_id: 0,
mean_m3s: 50.0,
std_m3s: 0.0,
ar_coefficients: vec![],
residual_std_ratio: 0.0,
},
InflowModel {
hydro_id: h2_id,
stage_id: 0,
mean_m3s: 50.0,
std_m3s: 0.0,
ar_coefficients: vec![],
residual_std_ratio: 0.0,
},
];
let load_models = vec![LoadModel {
bus_id: EntityId(1),
stage_id: 0,
mean_mw: 100.0,
std_mw: 0.0,
}];
let constraint = GenericConstraint {
id: EntityId(100),
name: "gc_sum_gen".to_string(),
description: None,
expression: ConstraintExpression {
terms: vec![
LinearTerm {
coefficient: 1.0,
variable: VariableRef::HydroGeneration {
hydro_id: h1_id,
block_id: None,
},
},
LinearTerm {
coefficient: 1.0,
variable: VariableRef::HydroGeneration {
hydro_id: h2_id,
block_id: None,
},
},
],
},
sense: ConstraintSense::LessEqual,
slack: SlackConfig {
enabled: false,
penalty: None,
},
};
let id_map: HashMap<i32, usize> = [(100_i32, 0)].into_iter().collect();
let rows = vec![(100_i32, 0_i32, None::<i32>, 200.0_f64)];
let generic_bounds = ResolvedGenericConstraintBounds::new(&id_map, rows.into_iter());
let resolved_bounds = ResolvedBounds::new(
2, 0, 0,
0,
0,
1, default_hydro_bounds(),
ThermalStageBounds {
min_generation_mw: 0.0,
max_generation_mw: 0.0,
},
LineStageBounds {
direct_mw: 0.0,
reverse_mw: 0.0,
},
PumpingStageBounds {
min_flow_m3s: 0.0,
max_flow_m3s: 0.0,
},
ContractStageBounds {
min_mw: 0.0,
max_mw: 0.0,
price_per_mwh: 0.0,
},
);
let penalties = ResolvedPenalties::new(
2, 1, 0, 0, 1, default_hydro_penalties(),
BusStagePenalties { excess_cost: 0.0 },
LineStagePenalties { exchange_cost: 0.0 },
NcsStagePenalties {
curtailment_cost: 0.0,
},
);
let system = SystemBuilder::new()
.buses(vec![bus])
.hydros(hydros)
.stages(vec![stage])
.inflow_models(inflow_models)
.load_models(load_models)
.bounds(resolved_bounds)
.penalties(penalties)
.generic_constraints(vec![constraint])
.resolved_generic_bounds(generic_bounds)
.build()
.expect("two_hydro_system: valid");
let t = &build_templates_for(&system)[0];
let generic_row = t.num_rows - 1;
let turbine_h1_col = 7_usize;
let turbine_h2_col = 8_usize;
let entries_h1 = csc_entries_at(t, turbine_h1_col, generic_row);
assert!(
!entries_h1.is_empty(),
"no CSC entry found for H1 turbine at generic constraint row"
);
let total_h1: f64 = entries_h1.iter().sum();
assert!(
(total_h1 - prod_h1).abs() < f64::EPSILON,
"expected coefficient {prod_h1} for H1, got {total_h1}"
);
let entries_h2 = csc_entries_at(t, turbine_h2_col, generic_row);
assert!(
!entries_h2.is_empty(),
"no CSC entry found for H2 turbine at generic constraint row"
);
let total_h2: f64 = entries_h2.iter().sum();
assert!(
(total_h2 - prod_h2).abs() < f64::EPSILON,
"expected coefficient {prod_h2} for H2, got {total_h2}"
);
}
#[allow(clippy::cast_possible_wrap)]
fn one_bus_one_thermal_system(
thermal_entity_id: EntityId,
constraints: Vec<cobre_core::GenericConstraint>,
bounds: cobre_core::ResolvedGenericConstraintBounds,
) -> cobre_core::System {
use chrono::NaiveDate;
use cobre_core::entities::thermal::{Thermal, ThermalCostSegment};
use cobre_core::scenario::LoadModel;
use cobre_core::temporal::{
Block, BlockMode, NoiseMethod, ScenarioSourceConfig, Stage, StageRiskConfig,
StageStateConfig,
};
let bus = Bus {
id: EntityId(1),
name: "B1".to_string(),
deficit_segments: vec![DeficitSegment {
depth_mw: None,
cost_per_mwh: 500.0,
}],
excess_cost: 0.0,
};
let thermal = Thermal {
id: thermal_entity_id,
name: "T1".to_string(),
bus_id: EntityId(1),
min_generation_mw: 0.0,
max_generation_mw: 100.0,
cost_segments: vec![ThermalCostSegment {
capacity_mw: 100.0,
cost_per_mwh: 50.0,
}],
gnl_config: None,
entry_stage_id: None,
exit_stage_id: None,
};
let stage = Stage {
index: 0,
id: 0_i32,
start_date: NaiveDate::from_ymd_opt(2024, 1, 1).unwrap(),
end_date: NaiveDate::from_ymd_opt(2024, 2, 1).unwrap(),
season_id: None,
blocks: vec![Block {
index: 0,
name: "BLK0".to_string(),
duration_hours: 744.0,
}],
block_mode: BlockMode::Parallel,
state_config: StageStateConfig {
storage: false,
inflow_lags: false,
},
risk_config: StageRiskConfig::Expectation,
scenario_config: ScenarioSourceConfig {
branching_factor: 1,
noise_method: NoiseMethod::Saa,
},
};
let load_models = vec![LoadModel {
bus_id: EntityId(1),
stage_id: 0,
mean_mw: 100.0,
std_mw: 0.0,
}];
let resolved_bounds = ResolvedBounds::new(
0,
1,
0,
0,
0,
1,
default_hydro_bounds(),
ThermalStageBounds {
min_generation_mw: 0.0,
max_generation_mw: 100.0,
},
LineStageBounds {
direct_mw: 0.0,
reverse_mw: 0.0,
},
PumpingStageBounds {
min_flow_m3s: 0.0,
max_flow_m3s: 0.0,
},
ContractStageBounds {
min_mw: 0.0,
max_mw: 0.0,
price_per_mwh: 0.0,
},
);
let penalties = ResolvedPenalties::new(
0,
1,
0,
0,
1,
default_hydro_penalties(),
BusStagePenalties { excess_cost: 0.0 },
LineStagePenalties { exchange_cost: 0.0 },
NcsStagePenalties {
curtailment_cost: 0.0,
},
);
SystemBuilder::new()
.buses(vec![bus])
.thermals(vec![thermal])
.stages(vec![stage])
.load_models(load_models)
.bounds(resolved_bounds)
.penalties(penalties)
.generic_constraints(constraints)
.resolved_generic_bounds(bounds)
.build()
.expect("one_bus_one_thermal_system: valid")
}
fn minimal_template(
num_rows: usize,
num_cols: usize,
col_starts: Vec<i32>,
row_indices: Vec<i32>,
values: Vec<f64>,
row_lower: Vec<f64>,
row_upper: Vec<f64>,
) -> StageTemplate {
let num_nz = values.len();
StageTemplate {
num_cols,
num_rows,
num_nz,
col_starts,
row_indices,
values,
col_lower: vec![0.0; num_cols],
col_upper: vec![f64::INFINITY; num_cols],
objective: vec![0.0; num_cols],
row_lower,
row_upper,
n_state: 0,
n_transfer: 0,
n_dual_relevant: 0,
n_hydro: 0,
max_par_order: 0,
col_scale: Vec::new(),
row_scale: Vec::new(),
}
}
#[test]
fn row_scale_identity_for_uniform_matrix() {
let col_starts = vec![0, 2, 4];
let row_indices = vec![0, 1, 0, 1];
let values = vec![1.0, 1.0, 1.0, 1.0];
let scale = super::compute_row_scale(2, 2, &col_starts, &row_indices, &values);
assert_eq!(scale.len(), 2);
assert!(
(scale[0] - 1.0).abs() < 1e-15,
"row 0 scale should be 1.0, got {}",
scale[0]
);
assert!(
(scale[1] - 1.0).abs() < 1e-15,
"row 1 scale should be 1.0, got {}",
scale[1]
);
}
#[test]
fn row_scale_geometric_mean() {
let col_starts = vec![0, 1, 3];
let row_indices = vec![0, 0, 1];
let values = vec![1.0, 100.0, 4.0];
let scale = super::compute_row_scale(2, 2, &col_starts, &row_indices, &values);
assert_eq!(scale.len(), 2);
let expected_row0 = 1.0_f64 / (1.0_f64 * 100.0_f64).sqrt(); let expected_row1 = 1.0_f64 / (4.0_f64 * 4.0_f64).sqrt(); assert!(
(scale[0] - expected_row0).abs() < 1e-14,
"row 0 scale: expected {expected_row0}, got {}",
scale[0]
);
assert!(
(scale[1] - expected_row1).abs() < 1e-14,
"row 1 scale: expected {expected_row1}, got {}",
scale[1]
);
}
#[test]
fn apply_row_scale_scales_values_and_bounds() {
let col_starts = vec![0_i32, 1, 3];
let row_indices = vec![0_i32, 0, 1];
let values = vec![1.0_f64, 100.0, 4.0];
let row_lower = vec![-5.0_f64, 7.0];
let row_upper = vec![f64::INFINITY, 7.0];
let mut tmpl =
minimal_template(2, 2, col_starts, row_indices, values, row_lower, row_upper);
let row_scale = vec![0.1_f64, 0.25];
super::apply_row_scale(&mut tmpl, &row_scale);
assert!((tmpl.values[0] - 0.1).abs() < 1e-15, "value[0] wrong");
assert!((tmpl.values[1] - 10.0).abs() < 1e-15, "value[1] wrong");
assert!((tmpl.values[2] - 1.0).abs() < 1e-15, "value[2] wrong");
assert!(
(tmpl.row_lower[0] - (-0.5)).abs() < 1e-15,
"row_lower[0] wrong"
);
assert!(
tmpl.row_upper[0].is_infinite() && tmpl.row_upper[0] > 0.0,
"row_upper[0] must remain +inf"
);
assert!(
(tmpl.row_lower[1] - 1.75).abs() < 1e-15,
"row_lower[1] wrong"
);
assert!(
(tmpl.row_upper[1] - 1.75).abs() < 1e-15,
"row_upper[1] wrong"
);
assert_eq!(tmpl.col_lower, vec![0.0; 2]);
assert!(tmpl.col_upper[0].is_infinite());
assert!(tmpl.col_upper[1].is_infinite());
assert_eq!(tmpl.objective, vec![0.0; 2]);
}
#[test]
fn row_scale_empty_row_gets_one() {
let col_starts = vec![0_i32, 1];
let row_indices = vec![1_i32];
let values = vec![8.0_f64];
let scale = super::compute_row_scale(3, 1, &col_starts, &row_indices, &values);
assert_eq!(scale.len(), 3);
assert!(
(scale[0] - 1.0).abs() < 1e-15,
"empty row 0 scale should be 1.0"
);
let expected = 1.0_f64 / 8.0;
assert!(
(scale[1] - expected).abs() < 1e-15,
"row 1 scale: expected {expected}, got {}",
scale[1]
);
assert!(
(scale[2] - 1.0).abs() < 1e-15,
"empty row 2 scale should be 1.0"
);
}
}