microgrid.m - master - Moustafa/Microgrid

microgrid.m
2 049 строк98 KB
Moustafa
create: Load_and_train_lstm.m, microgrid.m, microgrid_data.m, Microgrid_lstm_predict.m, Microgrid_lstm_train.m, P real data.docx, P real matrix.docx
25 апр 2026, 17:25
Верифицирован
25 апр 2026, 17:251d85e06
function microgrid()
% =========================================================================
%  microgrid — ПОЛНАЯ ВЕРСИЯ С ПОДДЕРЖКОЙ ВНЕШНЕГО УПРАВЛЕНИЯ
%  ВЕРСИЯ 22.2 — Все инженерные и экономические исправления:
%
%  [FIX-A] Верхний предел базовой нагрузки из db.Const (IEC 62898-1:2017)
%          Limits=[1, db.Const.microgrid_max_kw] вместо [1, Inf]
%  [FIX-B] update_caps_for_city(): реалистичный подбор N параллельных ДГУ,
%          ограничения ВЭС/СЭС по классам оборудования (ISO 8528-1:2018,
%          IEC 61400-1:2019)
%  [FIX-C] tiered_capex(): многоуровневый CAPEX по масштабу установки
%          (IRENA 2023, BloombergNEF 2024)
%  [FIX-D] Ставка дисконтирования 25% (ЦБ РФ 21% + 4% риск, апрель 2026)
%          — теперь только в microgrid_data.m (принцип DRY)
%  [FIX-E] VOLL (стоимость недопоставки): 350 руб/кВт·ч (ФАС РФ 2023)
%          учитывается в LCOE при наличии дефицита
%  [FIX-DRY] recalc_economics(): все константы из db.Const, не хардкод
%  [FIX-LSTM] LSTM — только для управления BESS (upd. bess_reserve),
%             P_load всегда = P_load_det (реальный профиль)
%  [FIX-FORMAT] format_power/format_energy: единицы автоматически кВт/МВт
%  [FIX-WARN] Предупреждение при нагрузке >microgrid_warn_kw
%  [FIX-1]  NPV: раздельная индексация топлива/O&M
%  [FIX-2]  BESS заряд ограничен p_bess_val (IEC 62619:2022)
%  [FIX-6]  Температурный коэффициент СЭС: -0.0040/°C (IEC 61215:2021)
%  [FIX-7]  bess_repl_year = 8 лет из db.Const (IEC 62619:2022)
%  [FIX-8]  Остаточная стоимость с линейной амортизацией (IRENA 2023)
%  [FIX-9]  Вейбулл параметр по регионам (Росгидромет 2022)
%  [FIX-12] Полный список праздников РФ (ТК РФ, ст.112) из db.Const
%
%  Источники:
%    IEC 62898-1:2017 — Microgrids, Part 1
%    IEC 62619:2022 — Lithium battery safety
%    ISO 8528-1:2018 — Generating sets
%    IEC 61400-1:2019 — Wind energy
%    IRENA "Renewable Power Generation Costs 2023"
%    BloombergNEF "New Energy Outlook 2024"
%    Банк России, ключевая ставка апрель 2026
%    ФАС РФ, VOLL, Методические рекомендации 2023
% =========================================================================

close all;
fig = uifigure('Name','Microgrid System v22.2','Position',[30 30 1600 950]);

% =========================================================================
% 1. БАЗА ДАННЫХ
% =========================================================================
db = microgrid_data();

% =========================================================================
% 2. СТРУКТУРА ДАННЫХ
% =========================================================================
setappdata(0, 'MicrogridFigure', fig);
modelData = struct('pv_opt',0,'wind_opt',0,'bess_p_opt',0,'bess_e_opt',0,...
    'dg1_opt',0,'dg2_opt',0,'dg_n_units',1,...
    'sources',struct('ses',true,'ves',true,'dgu1',true,'dgu2',true,'bess',true));
setappdata(fig,'modelData',modelData);
setappdata(fig,'dg_hours',[0,0]);

% =========================================================================
% 3. ИНТЕРФЕЙС
% =========================================================================
uilabel(fig,'Position',[15 895 200 20],'Text','Регион:','FontWeight','bold');
regDrop = uidropdown(fig,'Position',[15 872 210 22],...
    'Items',db.RegionNames,'ItemsData',db.RegionKeys);

uilabel(fig,'Position',[15 845 200 20],'Text','Город:','FontWeight','bold');
cityDrop = uidropdown(fig,'Position',[15 822 210 22],...
    'Items',db.Cities.(db.RegionKeys{1}));

uilabel(fig,'Position',[15 795 200 20],'Text','Базовая нагрузка (кВт):','FontWeight','bold');
% [FIX-A] Верхний предел из db.Const.microgrid_max_kw (IEC 62898-1:2017)
% Не Inf — только реальный диапазон микросети
loadField = uieditfield(fig,'numeric','Position',[15 772 210 22],...
    'Value',500,...
    'Limits',[db.Const.microgrid_min_kw, db.Const.microgrid_max_kw]);

uilabel(fig,'Position',[15 745 200 20],'Text','Дата:','FontWeight','bold');
datePick = uidatepicker(fig,'Position',[15 722 210 22],'Value',datetime('today'));

uilabel(fig,'Position',[15 695 200 20],'Text','Тип нагрузки:','FontWeight','bold');
loadTypeDrop = uidropdown(fig,'Position',[15 672 210 22],...
    'Items',{'Смешанный','Промышленный','Жилой'},...
    'ItemsData',{'mixed','industrial','residential'});

% --- Метки мощностей ---
uilabel(fig,'Position',[15 477 40 18],'Text','СЭС:','FontSize',9);
pvLabel = uilabel(fig,'Position',[55 477 72 18],'Text','0 кВт',...
    'FontSize',9,'FontWeight','bold','Tag','pvLabel');
setappdata(fig,'pvLabel',pvLabel);

uilabel(fig,'Position',[132 477 35 18],'Text','ВЭС:','FontSize',9);
windLabel = uilabel(fig,'Position',[168 477 67 18],'Text','0 кВт',...
    'FontSize',9,'FontWeight','bold','Tag','windLabel');
setappdata(fig,'windLabel',windLabel);

uilabel(fig,'Position',[15 455 55 18],'Text','BESS кВт:','FontSize',9);
bessPLabel = uilabel(fig,'Position',[70 455 45 18],'Text','0 кВт',...
    'FontSize',9,'FontWeight','bold','Tag','bessPLabel');
setappdata(fig,'bessPLabel',bessPLabel);

uilabel(fig,'Position',[118 455 60 18],'Text','BESS кВт·ч:','FontSize',9);
bessELabel = uilabel(fig,'Position',[178 455 53 18],'Text','0 кВт·ч',...
    'FontSize',9,'FontWeight','bold','Tag','bessELabel');
setappdata(fig,'bessELabel',bessELabel);

uilabel(fig,'Position',[15 433 55 18],'Text','ДГУ1 кВт:','FontSize',9);
dg1Label = uilabel(fig,'Position',[70 433 60 18],'Text','0 кВт',...
    'FontSize',9,'FontWeight','bold','Tag','dg1Label');
setappdata(fig,'dg1Label',dg1Label);

uilabel(fig,'Position',[135 433 55 18],'Text','ДГУ2 кВт:','FontSize',9);
dg2Label = uilabel(fig,'Position',[190 433 30 18],'Text','0',...
    'FontSize',9,'FontWeight','bold','Tag','dg2Label');
setappdata(fig,'dg2Label',dg2Label);

% --- Кнопки источников ---
btnPV = uibutton(fig,'state','Position',[15 573 100 22],'Text','СЭС ON',...
    'Value',true,'BackgroundColor',[0.2 0.75 0.2],'FontColor','w','FontWeight','bold','Tag','btnPV');
setappdata(fig,'btnPV',btnPV);

btnWind = uibutton(fig,'state','Position',[123 573 100 22],'Text','ВЭС ON',...
    'Value',true,'BackgroundColor',[0.2 0.75 0.2],'FontColor','w','FontWeight','bold','Tag','btnWind');
setappdata(fig,'btnWind',btnWind);

btnDG1 = uibutton(fig,'state','Position',[15 548 100 22],'Text','ДГУ1 ON',...
    'Value',true,'BackgroundColor',[0.2 0.75 0.2],'FontColor','w','FontWeight','bold','Tag','btnDG1');
setappdata(fig,'btnDG1',btnDG1);

btnDG2 = uibutton(fig,'state','Position',[123 548 100 22],'Text','ДГУ2 ON',...
    'Value',true,'BackgroundColor',[0.2 0.75 0.2],'FontColor','w','FontWeight','bold','Tag','btnDG2');
setappdata(fig,'btnDG2',btnDG2);

btnBESS = uibutton(fig,'state','Position',[15 522 210 22],'Text','BESS ON',...
    'Value',true,'BackgroundColor',[0.2 0.75 0.2],'FontColor','w','FontWeight','bold','Tag','btnBESS');
setappdata(fig,'btnBESS',btnBESS);

runBtn = uibutton(fig,'Position',[15 622 210 46],'Text','▶  МОДЕЛИРОВАТЬ',...
    'BackgroundColor',[0.1 0.4 0.7],'FontColor','w','FontWeight','bold','FontSize',13,'Tag','runBtn');
setappdata(fig,'runBtn',runBtn);

optBtn = uibutton(fig,'Position',[15 405 210 24],'Text','🔧  Авто-оптимизация ВИЭ',...
    'BackgroundColor',[0.45 0.15 0.55],'FontColor','w','FontWeight','bold','FontSize',10,'Tag','optBtn');
setappdata(fig,'optBtn',optBtn);

resetSOCBtn = uibutton(fig,'Position',[15 378 210 22],'Text','⟳ Сбросить SOC (80%)',...
    'BackgroundColor',[0.3 0.5 0.7],'FontColor','w','FontWeight','bold','FontSize',10,'Tag','resetSOCBtn');

reportBtn = uibutton(fig,'Position',[15 162 210 54],...
    'Text',sprintf('📊 Сравнит.\nОтчёт (экономика)'),...
    'BackgroundColor',[0.85 0.45 0.05],'FontColor','w','FontWeight','bold','FontSize',11,'Tag','reportBtn');
setappdata(fig,'reportBtn',reportBtn);

dguReportBtn = uibutton(fig,'Position',[15 108 210 54],...
    'Text',sprintf('📊 Отчёт по\nДГУ1/ДГУ2'),...
    'BackgroundColor',[0.2 0.5 0.7],'FontColor','w','FontWeight','bold','FontSize',11,'Tag','dguReportBtn');
setappdata(fig,'dguReportBtn',dguReportBtn);
dguReportBtn.ButtonPushedFcn = @(btn,e) open_dgu_report();

dguCharsBtn = uibutton(fig,'Position',[15 54 210 54],...
    'Text',sprintf('📈 Характеристики\nДГУ (B(P) и α(P))'),...
    'BackgroundColor',[0.3 0.6 0.4],'FontColor','w','FontWeight','bold','FontSize',11,'Tag','dguCharsBtn');
dguCharsBtn.ButtonPushedFcn = @(btn,e) open_dgu_chars();

diagArea = uitextarea(fig,'Position',[15 222 210 157],...
    'Editable',false,'FontSize',9,'BackgroundColor',[0.97 0.97 0.97],'Tag','diagArea');
setappdata(fig,'diagArea',diagArea);

uilabel(fig,'Position',[240 8 200 18],'Text','Microgrid AI v22.2',...
    'FontSize',7,'FontColor',[0.5 0.5 0.5]);

% =========================================================================
% 4. ФУНКЦИИ ВНЕШНЕГО УПРАВЛЕНИЯ
% =========================================================================
function success = setSourcePower(source, value)
    success = false;
    if ~isnumeric(value) || ~isfinite(value) || value < 0, return; end
    try
        md = getappdata(fig,'modelData');
        switch source
            case 'ses',    md.pv_opt = value;     pvLabel.Text  = format_power(value);
            case 'ves',    md.wind_opt = value;    windLabel.Text = format_power(value);
            case 'dgu1',   md.dg1_opt = value;     dg1Label.Text = format_power(value);
            case 'dgu2',   md.dg2_opt = value;     dg2Label.Text = format_power(value);
            case 'bess',   md.bess_p_opt = value;  bessPLabel.Text = format_power(value);
            case 'bess_e', md.bess_e_opt = value;  bessELabel.Text = format_energy(value);
            otherwise, return;
        end
        setappdata(fig,'modelData',md); run_logic(); success = true;
    catch e
        disp(['Ошибка setSourcePower: ' e.message]);
    end
end

function success = setSourceState(source, state)
    success = false;
    try
        md = getappdata(fig,'modelData');
        switch source
            case 'ses',  btnPV.Value=state;   btnPV.Text=iff(state,'СЭС ON','СЭС OFF');   btnPV.BackgroundColor=iff(state,[0.2 0.75 0.2],[0.8 0.2 0.2]);   md.sources.ses=state;
            case 'ves',  btnWind.Value=state;  btnWind.Text=iff(state,'ВЭС ON','ВЭС OFF'); btnWind.BackgroundColor=iff(state,[0.2 0.75 0.2],[0.8 0.2 0.2]); md.sources.ves=state;
            case 'dgu1', btnDG1.Value=state;   btnDG1.Text=iff(state,'ДГУ1 ON','ДГУ1 OFF'); btnDG1.BackgroundColor=iff(state,[0.2 0.75 0.2],[0.8 0.2 0.2]); md.sources.dgu1=state;
            case 'dgu2', btnDG2.Value=state;   btnDG2.Text=iff(state,'ДГУ2 ON','ДГУ2 OFF'); btnDG2.BackgroundColor=iff(state,[0.2 0.75 0.2],[0.8 0.2 0.2]); md.sources.dgu2=state;
            case 'bess', btnBESS.Value=state;  btnBESS.Text=iff(state,'BESS ON','BESS OFF'); btnBESS.BackgroundColor=iff(state,[0.2 0.75 0.2],[0.8 0.2 0.2]); md.sources.bess=state;
            otherwise, return;
        end
        setappdata(fig,'modelData',md); run_logic(); success = true;
    catch e
        disp(['Ошибка setSourceState: ' e.message]);
    end
end

function success = runExternal()
    success = false;
    try, run_logic(); success = true;
    catch e, disp(['Ошибка runExternal: ' e.message]); end
end

function [powers, sources] = getModelData()
    md = getappdata(fig,'modelData');
    powers = struct('ses',md.pv_opt,'ves',md.wind_opt,'bess',md.bess_p_opt,...
        'bess_e',md.bess_e_opt,'dgu1',md.dg1_opt,'dgu2',md.dg2_opt);
    sources = struct('ses',btnPV.Value,'ves',btnWind.Value,'dgu1',btnDG1.Value,...
        'dgu2',btnDG2.Value,'bess',btnBESS.Value);
end

setappdata(fig,'setSourcePowerFcn',@setSourcePower);
setappdata(fig,'setSourceStateFcn',@setSourceState);
setappdata(fig,'runExternalFcn',@runExternal);
setappdata(fig,'getModelDataFcn',@getModelData);

% =========================================================================
% 5. ОСИ ГРАФИКОВ
% =========================================================================
ax1 = uiaxes(fig,'Position',[240 720 310 200]);
ax2 = uiaxes(fig,'Position',[565 720 310 200]);
ax4 = uiaxes(fig,'Position',[890 720 310 200]);
ax5 = uiaxes(fig,'Position',[1196 720 379 200]);
ax3 = uiaxes(fig,'Position',[240 30 1335 660]);

fig.WindowButtonDownFcn   = @on_drag_start;
fig.WindowButtonMotionFcn = @on_drag_move;
fig.WindowButtonUpFcn     = @on_drag_end;

% =========================================================================
% 6. ОБРАБОТЧИКИ
% =========================================================================
btnPV.ValueChangedFcn   = @(b,e) toggle_source(b,'СЭС');
btnWind.ValueChangedFcn = @(b,e) toggle_source(b,'ВЭС');
btnDG1.ValueChangedFcn  = @(b,e) toggle_source(b,'ДГУ1');
btnDG2.ValueChangedFcn  = @(b,e) toggle_source(b,'ДГУ2');
btnBESS.ValueChangedFcn = @(b,e) toggle_source(b,'BESS');

regDrop.ValueChangedFcn      = @(dd,e) updateCities();
cityDrop.ValueChangedFcn     = @(dd,e) on_city_change();
datePick.ValueChangedFcn     = @(dd,e) run_logic();
loadField.ValueChangedFcn    = @(dd,e) on_load_changed();
loadTypeDrop.ValueChangedFcn = @(dd,e) run_logic();
runBtn.ButtonPushedFcn       = @(btn,e) run_simulate_all_days();
reportBtn.ButtonPushedFcn    = @(btn,e) open_report();
optBtn.ButtonPushedFcn       = @(btn,e) run_opt_then_simulate();
resetSOCBtn.ButtonPushedFcn  = @(btn,e) reset_soc();

% Переменные
lastShownDeficit = 0;
lastResults = [];
allDaysData = struct('P_load',[],'P_PV',[],'P_Wind',[],'P_DG1',[],'P_DG2',[],...
    'P_BESS',[],'SOC',[],'P_deficit',[],'P_surplus',[],'dates',{{}});
dragState = struct('active',false,'startPixX',0,'startXLim',[-0.5 23.5],'dir',0);
batch_mode = false;

% ИНИЦИАЛИЗАЦИЯ
md = getappdata(fig,'modelData');
md.pv_opt=100; md.wind_opt=45; md.bess_p_opt=80;
md.bess_e_opt=400; md.dg1_opt=325; md.dg2_opt=275; md.dg_n_units=1;
setappdata(fig,'modelData',md);
pvLabel.Text    = format_power(md.pv_opt);
windLabel.Text  = format_power(md.wind_opt);
bessPLabel.Text = format_power(md.bess_p_opt);
bessELabel.Text = format_energy(md.bess_e_opt);
dg1Label.Text   = format_power(md.dg1_opt);
dg2Label.Text   = format_power(md.dg2_opt);

% =========================================================================
% 7. ВСПОМОГАТЕЛЬНЫЕ ФУНКЦИИ
% =========================================================================

% [FIX-FORMAT] Автоматические единицы кВт/МВт
function s = format_power(kw)
    if kw >= 1000
        s = sprintf('%.2f МВт', kw/1000);
    else
        s = sprintf('%d кВт', round(kw));
    end
end

function s = format_energy(kwh)
    if kwh >= 1000
        s = sprintf('%.2f МВт·ч', kwh/1000);
    else
        s = sprintf('%d кВт·ч', round(kwh));
    end
end

% [FIX-C] Многоуровневый CAPEX по масштабу установки
% Источник: IRENA "Renewable Power Generation Costs 2023";
%           BloombergNEF "New Energy Outlook 2024"
function cost = tiered_capex(capacity, tiers)
    % tiers: матрица [до_единицы, руб/единица]
    cost = 0;
    if capacity <= 0, return; end
    for i = 1:size(tiers,1)
        if capacity <= tiers(i,1)
            cost = capacity * tiers(i,2);
            return;
        end
    end
    cost = capacity * tiers(end,2);
end

function updateCities()
    cityDrop.Items = db.Cities.(regDrop.Value);
    cityDrop.Value = cityDrop.Items{1};
    update_caps_for_city();
    run_logic();
end

function on_city_change()
    update_caps_for_city();
    run_logic();
end

% [FIX-A] [FIX-WARN] Валидация нагрузки при изменении поля
function on_load_changed()
    v = loadField.Value;
    if v > db.Const.microgrid_warn_kw
        lines = {};
        lines{1} = sprintf('⚠ НАГРУЗКА %.0f кВт', v);
        lines{2} = sprintf('> предела микросети (%d кВт)', db.Const.microgrid_warn_kw);
        lines{3} = '→ Рассчитывается как параллельные';
        lines{4} = 'агрегаты (N×единичная мощность).';
        lines{5} = 'Физически реализуемо.';
        diagArea.Value = lines;
    end
    update_caps_for_city();
    run_logic();
end

function reset_soc()
    setappdata(fig,'last_soc_value',80);
    run_logic();
end

% [FIX-B] Реалистичный подбор оборудования с учётом классов и ограничений
% Источники:
%   ISO 8528-1:2018 — Rating and performance of generating sets
%   IEC 61400-1:2019 — Design requirements for wind turbines
%   Методические рекомендации Минэнерго РФ 2022
%   EPRI 3002019190 (2020) — оптимальная BESS = 20-25% P_base
function update_caps_for_city()
    md = getappdata(fig,'modelData');
    base_load = max(loadField.Value, 1);
    vp = get_vie_profile(regDrop.Value, cityDrop.Value);

    % ── ДГУ: параллельная работа N агрегатов ─────────────────────────
    % ISO 8528-1:2018 — max единичная мощность серийного ДГУ ≈ 2500 кВт
    dg_single_max = db.Const.dg_max_single_kw;   % 2500 кВт
    dg_total_need_1 = round(base_load * 0.65);
    dg_total_need_2 = round(base_load * 0.55);

    % Количество агрегатов каждого типа
    n1 = max(1, ceil(dg_total_need_1 / dg_single_max));
    n1 = min(n1, db.Const.dg_max_units);  % не более 4 агрегатов
    n2 = max(1, ceil(dg_total_need_2 / dg_single_max));
    n2 = min(n2, db.Const.dg_max_units);

    % Единичная мощность: округляем до 50 кВт
    dg1_single = min(round(dg_total_need_1 / n1 / 50) * 50, dg_single_max);
    dg1_single = max(dg1_single, db.Const.dg_min_single_kw);
    dg2_single = min(round(dg_total_need_2 / n2 / 50) * 50, dg_single_max);
    dg2_single = max(dg2_single, db.Const.dg_min_single_kw);

    md.dg1_opt = dg1_single;
    md.dg2_opt = dg2_single;
    md.dg_n_units = max(n1, n2);   % сохраняем для отображения

    % ── СЭС: практический предел площадки ──────────────────────────────
    pv_rec = max(round(base_load * vp.pv_mult_rec / 50) * 50, 0);
    pv_max = min(base_load * 2.0, db.Const.pv_max_kw);
    md.pv_opt = min(pv_rec, pv_max);

    % ── ВЭС: класс малых ВЭУ (IEC 61400-1:2019) ────────────────────────
    wind_rec = max(round(base_load * vp.wind_mult_rec / 10) * 10, 0);
    wind_max = min(base_load * 1.5, db.Const.wind_max_total_kw);
    md.wind_opt = min(wind_rec, wind_max);

    % ── BESS: EPRI 3002019190 (2020) оптимум 20-25% ────────────────────
    md.bess_p_opt = max(round(base_load * 0.22 / 10) * 10, 10);
    md.bess_p_opt = min(md.bess_p_opt, db.Const.bess_max_kw);
    md.bess_e_opt = max(round(md.bess_p_opt * 5 / 50) * 50, 50);
    md.bess_e_opt = min(md.bess_e_opt, db.Const.bess_max_kwh);

    setappdata(fig,'modelData',md);

    % Обновление меток с единицами
    pvLabel.Text    = format_power(md.pv_opt);
    windLabel.Text  = format_power(md.wind_opt);
    bessPLabel.Text = format_power(md.bess_p_opt);
    bessELabel.Text = format_energy(md.bess_e_opt);
    % Для ДГУ показываем единичную мощность × количество агрегатов
    if md.dg_n_units > 1
        dg1Label.Text = sprintf('%s×%d', format_power(md.dg1_opt), n1);
        dg2Label.Text = sprintf('%s×%d', format_power(md.dg2_opt), n2);
    else
        dg1Label.Text = format_power(md.dg1_opt);
        dg2Label.Text = format_power(md.dg2_opt);
    end
end

function run_simulate_all_days()
    original_date = datePick.Value;
    orig_dn = datenum(original_date);
    N_future = 7;
    all_dates_to_calc = {};
    for di = 1:length(allDaysData.dates)
        all_dates_to_calc{end+1} = allDaysData.dates{di};
    end
    if isempty(all_dates_to_calc) || ...
       ~any(cellfun(@(d) abs(datenum(d)-orig_dn)<0.5, all_dates_to_calc))
        all_dates_to_calc{end+1} = original_date;
    end
    for fi = 1:N_future
        fd = original_date + caldays(fi);
        fd_dn = datenum(fd);
        if ~any(cellfun(@(d) abs(datenum(d)-fd_dn)<0.5, all_dates_to_calc))
            all_dates_to_calc{end+1} = fd;
        end
    end
    dns = cellfun(@(d) datenum(d), all_dates_to_calc);
    [~,si] = sort(dns);
    all_dates_to_calc = all_dates_to_calc(si);
    allDaysData = struct('P_load',[],'P_PV',[],'P_Wind',[],'P_DG1',[],'P_DG2',[],...
        'P_BESS',[],'SOC',[],'P_deficit',[],'P_surplus',[],'dates',{{}});
    setappdata(fig,'last_soc_value',80);
    setappdata(fig,'dg_hours',[0,0]);
    setappdata(fig,'bess_soh',1.0);
    datePick.ValueChangedFcn = [];
    batch_mode = true;
    dayResultsCache = {};
    datePick.Visible = 'off';
    try
        for di = 1:length(all_dates_to_calc)
            datePick.Value = all_dates_to_calc{di};
            bess_deg_per_day = db.Const.bess_deg_rate / 365;
            current_soh = max(0.80, 1.0 - bess_deg_per_day * (di-1));
            setappdata(fig,'bess_soh',current_soh);
            run_logic();
            dayResultsCache{end+1} = lastResults;
        end
    catch simErr
        batch_mode = false;
        datePick.Value = original_date;
        datePick.Visible = 'on';
        datePick.ValueChangedFcn = @(dd,e) run_logic();
        rethrow(simErr);
    end
    batch_mode = false;
    datePick.Value = original_date;
    datePick.Visible = 'on';
    datePick.ValueChangedFcn = @(dd,e) run_logic();
    for di = 1:length(all_dates_to_calc)
        if abs(datenum(all_dates_to_calc{di})-orig_dn) < 0.5
            if di <= length(dayResultsCache) && ~isempty(dayResultsCache{di})
                lastResults = dayResultsCache{di};
                recalc_economics();
            end
            break;
        end
    end
    update_graphs_with_new_data();
    drawnow;
end

function run_opt_then_simulate()
    run_optimize();
    run_simulate_all_days();
end

function season = get_season(dateVal)
    m = month(dateVal);
    if m>=12||m<=2, season='Winter';
    elseif m<=5, season='Spring';
    elseif m<=8, season='Summer';
    else, season='Autumn';
    end
end

function dtype = get_day_type(dateVal)
    dow = weekday(dateVal);
    if dow==1||dow==7, dtype='Weekend'; else, dtype='Weekday'; end
end

function key = city_key(cityName)
    if isKey(db.CityKey,cityName), key=db.CityKey(cityName);
    else, key='Moscow'; end
end

function k = holiday_multiplier(city, dateVal)
    % [FIX-12] Праздники из db.Const.holidays (ТК РФ ст.112)
    k=1.0;
    dayNum=day(dateVal); monthNum=month(dateVal);
    ckey=city_key(city);
    if isfield(db.CityType,ckey), ctype=db.CityType.(ckey); else, ctype='mixed'; end
    hb=0;
    htbl = db.Const.holidays;
    for hi=1:size(htbl,1)
        if monthNum==htbl(hi,1) && dayNum==htbl(hi,2)
            hb=htbl(hi,3); break;
        end
    end
    if hb > 0
        switch ctype
            case 'industrial',  k = 1.0 - hb * 0.6;
            case 'residential', k = 1.0 + hb;
            case 'mixed',       k = 1.0 + hb * 0.3;
        end
    end
end

function lat = get_latitude_safe(cityName)
    ckey = city_key(cityName);
    if isfield(db.Latitude,ckey), lat=db.Latitude.(ckey);
    else, lat=55.0; end
end

function P_PV_out = calc_pv(t, p_pv_val, city, dateVal, varargin)
    lat = get_latitude_safe(city);
    doy=day(dateVal,'dayofyear');
    decl=23.45*sin(deg2rad(360/365*(doy-81)));
    cos_omega=max(-1,min(1,-tan(deg2rad(lat))*tan(deg2rad(decl))));
    daylight_h=(2/15)*rad2deg(acos(cos_omega));
    t_rise=12-daylight_h/2;
    if nargin>=5, ckey_reg=varargin{1}; else, ckey_reg=''; end
    cur_season=get_season(dateVal);
    if isfield(db.Econ,ckey_reg)
        re_pv=db.Econ.(ckey_reg);
        switch cur_season
            case 'Winter', solar_frac_k=re_pv.solar_frac_Win; T_air_base=re_pv.T_Win;
            case 'Spring', solar_frac_k=re_pv.solar_frac_Spr; T_air_base=re_pv.T_Spr;
            case 'Summer', solar_frac_k=re_pv.solar_frac_Sum; T_air_base=re_pv.T_Sum;
            case 'Autumn', solar_frac_k=re_pv.solar_frac_Aut; T_air_base=re_pv.T_Aut;
            otherwise,     solar_frac_k=0.65;            T_air_base=10;
        end
    else
        switch cur_season
            case 'Summer', solar_frac_k=0.85; T_air_base=20;
            case 'Winter', solar_frac_k=0.55; T_air_base=-8;
            otherwise,     solar_frac_k=0.70; T_air_base=8;
        end
    end
    pv_system_losses=db.Const.pv_system_losses;
    diffuse_fraction=0.20*(1-solar_frac_k);
    P_PV_out=zeros(1,24);
    for h=1:24
        th=t(h)+0.5;
        if th>=t_rise && th<=(t_rise+daylight_h)
            angle=pi*(th-t_rise)/daylight_h;
            irr_direct=sin(angle);
            irr_total=min(irr_direct*solar_frac_k + irr_direct*diffuse_fraction, 1.0);
            NOCT_delta=(45-20)/800*1000;
            T_cell_h=T_air_base + NOCT_delta*irr_direct;
            % [FIX-6] IEC 61215:2021 — монокристаллический Si: -0.40%/°C
            k_temp_h=max(0, 1 - db.Const.pv_temp_coeff*(T_cell_h-25));
            P_PV_out(h)=p_pv_val*irr_total*k_temp_h*pv_system_losses;
        end
    end
    P_PV_out=max(P_PV_out.*(1+0.04*randn(1,24)),0);
end

function P_Wind_out = calc_wind(p_wind_val, season, t_arr, varargin)
    v_cutin=3; v_rated=12; v_cutout=25;
    hub_k=(80/10)^0.14;
    % [FIX-9] Параметр Вейбулла по классу региона (Росгидромет 2022)
    weibull_k=db.Const.weibull_k_default;
    if nargin>=4 && isfield(db.WeibullK,varargin{1})
        k_class=db.WeibullK.(varargin{1});
        weibull_k=db.WeibullKValue.(k_class);
    end
    if nargin>=4 && isfield(db.Econ,varargin{1})
        re_w=db.Econ.(varargin{1});
        switch season
            case 'Winter', v_mean=re_w.v_Win*hub_k;
            case 'Spring', v_mean=re_w.v_Spr*hub_k;
            case 'Summer', v_mean=re_w.v_Sum*hub_k;
            case 'Autumn', v_mean=re_w.v_Aut*hub_k;
            otherwise,     v_mean=re_w.v_Sum*hub_k;
        end
    else
        switch season
            case 'Winter', v_mean=8.5*hub_k; case 'Summer', v_mean=5.5*hub_k;
            otherwise,     v_mean=7.0*hub_k;
        end
    end
    wind_elec_losses=0.94;
    diurnal=1-0.25*sin(pi*(t_arr+6)/24);
    noise=ar1_noise(24,0.15,0.85);
    P_Wind_out=zeros(1,24);
    denom3=v_rated^3-v_cutin^3;
    for h=1:24
        u_rand=max(rand(),1e-6);
        v_rayleigh_factor=(-log(max(u_rand,1e-9)))^(1/weibull_k);
        v=max(0, v_mean*diurnal(h)*v_rayleigh_factor+noise(h));
        if v<v_cutin || v>v_cutout
            P_Wind_out(h)=0;
        elseif v<v_rated
            P_Wind_out(h)=p_wind_val*(v^3-v_cutin^3)/denom3*wind_elec_losses;
        else
            P_Wind_out(h)=p_wind_val*wind_elec_losses;
        end
    end
    P_Wind_out=min(P_Wind_out, p_wind_val);
end

function noise = ar1_noise(n, sigma, phi)
    noise=zeros(1,n); noise(1)=sigma*randn();
    for i=2:n, noise(i)=phi*noise(i-1)+sigma*randn(); end
    actual_std=std(noise)+1e-9;
    % [FIX шум] Центрируем после нормализации для устранения смещения
    noise=noise/actual_std*sigma;
    noise=noise-mean(noise);  % убираем систематическое смещение
end

function toggle_source(btn, name)
    if btn.Value
        btn.Text=[name ' ON']; btn.BackgroundColor=[0.2 0.75 0.2];
    else
        btn.Text=[name ' OFF']; btn.BackgroundColor=[0.8 0.2 0.2];
    end
    md=getappdata(fig,'modelData');
    switch name
        case 'СЭС',  md.sources.ses=btn.Value;
        case 'ВЭС',  md.sources.ves=btn.Value;
        case 'ДГУ1', md.sources.dgu1=btn.Value;
        case 'ДГУ2', md.sources.dgu2=btn.Value;
        case 'BESS', md.sources.bess=btn.Value;
    end
    setappdata(fig,'modelData',md);
    run_logic();
end

function price = get_city_fuel(cityName, regKey)
    ckey=city_key(cityName);
    if isKey(db.FuelCity,ckey), price=db.FuelCity(ckey);
    elseif isfield(db.Econ,regKey), price=db.Econ.(regKey).fuel;
    else, price=78; end
end

function vp = get_vie_profile(regKey, cityName)
    hub_k=(80/10)^0.14;
    if isfield(db.Econ,regKey)
        re=db.Econ.(regKey);
        cloud_avg=(re.solar_frac_Win+re.solar_frac_Spr+re.solar_frac_Sum+re.solar_frac_Aut)/4;
        v_avg_10m=(re.v_Win+re.v_Spr+re.v_Sum+re.v_Aut)/4;
    else
        cloud_avg=0.60; v_avg_10m=5.0;
    end
    v_avg_hub=v_avg_10m*hub_k;
    lat=get_latitude_safe(cityName);
    lat_pv=max(0.05, 1.0-max(lat-35,0)/40);
    pv_score=cloud_avg*lat_pv*100;
    if v_avg_hub>=12, wind_score=100;
    elseif v_avg_hub>3
        wind_score=(v_avg_hub^3-27)/(1728-27)*100;
        wind_score=max(0,min(100,wind_score));
    else, wind_score=0; end
    total=pv_score+wind_score+1e-9;
    pv_frac=pv_score/total; wind_frac=wind_score/total;
    pv_viable=(pv_score>=20);
    wind_viable=(v_avg_hub>=5.0);
    if pv_viable && ~wind_viable
        label=sprintf('☀ Только СЭС (ветер %.1f м/с)',v_avg_hub); category='pv';
    elseif wind_viable && ~pv_viable
        label=sprintf('💨 Только ВЭС (солнце %.0f%%)',pv_score); category='wind';
    elseif ~pv_viable && ~wind_viable
        label=sprintf('⚠ Нет ВИЭ (ветер %.1f м/с)',v_avg_hub); category='none';
    elseif pv_frac>=0.65
        label=sprintf('☀ СЭС доминирует %.0f%%',pv_frac*100); category='pv';
    elseif wind_frac>=0.65
        label=sprintf('💨 ВЭС доминирует %.0f%%',wind_frac*100); category='wind';
    else
        label='⚖ Баланс СЭС/ВЭС'; category='balanced';
    end
    if pv_viable && wind_viable
        pv_m=min(0.8+pv_frac*1.7, 2.5); wind_m=min(0.5+wind_frac*2.0, 2.5);
    elseif pv_viable
        pv_m=min(2.0+pv_frac*0.5, 2.5); wind_m=0;
    elseif wind_viable
        pv_m=0; wind_m=min(1.5+wind_frac*0.5, 2.5);
    else
        pv_m=0; wind_m=0.5;
    end
    vp=struct('pv_score',pv_score,'wind_score',wind_score,'pv_frac',pv_frac,...
        'wind_frac',wind_frac,'category',category,'label',label,'lat',lat,...
        'cloud_avg',cloud_avg,'v_hub',v_avg_hub,'pv_mult_rec',pv_m,...
        'wind_mult_rec',wind_m,'pv_viable',pv_viable,'wind_viable',wind_viable);
end

function run_optimize()
    optBtn.Text='⏳ Оптимизация...'; optBtn.BackgroundColor=[0.7 0.5 0.1]; drawnow;
    base_opt=max(loadField.Value,1);
    ckey_reg_opt=regDrop.Value;
    rng(42);
    seasons={'Winter','Spring','Summer','Autumn'};
    days_s=[91,92,92,90];
    ref_dates={datetime(2025,1,15),datetime(2025,4,15),datetime(2025,7,15),datetime(2025,10,15)};
    ltm=db.LoadType.(loadTypeDrop.Value);
    P_load_s=cell(1,4); P_PV_s=cell(1,4); P_Wind_s=cell(1,4);
    for si=1:4
        sn=seasons{si};
        bp=db.Profiles.(sn).Weekday.*ltm;
        bpm=max(bp); if bpm>0, bp=bp/bpm; end
        ln=1+ar1_noise(24,0.06,0.7);
        P_load_s{si}=max(base_opt*bp.*ln,0);
        P_PV_s{si}=calc_pv(0:23,1.0,cityDrop.Value,ref_dates{si},ckey_reg_opt);
        P_Wind_s{si}=calc_wind(1.0,sn,0:23,ckey_reg_opt);
    end
    vp=get_vie_profile(ckey_reg_opt,cityDrop.Value);
    pv_mult=[0.5,0.7,0.9,1.0,1.2,1.5,1.8,2.0];
    wind_mult=[0.3,0.5,0.7,0.9,1.0,1.2,1.5,1.8];
    bess_mult=[0.2,0.3,0.4,0.5,0.6,0.8,1.0,1.2];
    bess_h=[2,3,4,5,6,8];
    if ~vp.pv_viable,   pv_mult=[]; end
    if ~vp.wind_viable, wind_mult=[]; end

    % [FIX-B] Реалистичные ДГУ для оптимизации
    dg_single_max=db.Const.dg_max_single_kw;
    dg_total_1=round(base_opt*0.65);
    dg_total_2=round(base_opt*0.55);
    n1=min(max(1,ceil(dg_total_1/dg_single_max)),db.Const.dg_max_units);
    n2=min(max(1,ceil(dg_total_2/dg_single_max)),db.Const.dg_max_units);
    dg1_opt_val=min(round(dg_total_1/n1/50)*50, dg_single_max);
    dg2_opt_val=min(round(dg_total_2/n2/50)*50, dg_single_max);

    if isfield(db.Econ,ckey_reg_opt)
        fp_o=get_city_fuel(cityDrop.Value,ckey_reg_opt);
        et_o=db.Econ.(ckey_reg_opt).tariff;
    else
        fp_o=get_city_fuel(cityDrop.Value,ckey_reg_opt); et_o=5.5;
    end
    best_payback=Inf;
    best_pv=0; best_wind=0; best_bess_p=0; best_bess_e=0;
    for ipv=1:max(numel(pv_mult),1)
        pv_kw=0; if ~isempty(pv_mult), pv_kw=round(base_opt*pv_mult(ipv)); end
        pv_kw=min(pv_kw, db.Const.pv_max_kw);
        for iw=1:max(numel(wind_mult),1)
            wind_kw=0; if ~isempty(wind_mult), wind_kw=round(base_opt*wind_mult(iw)); end
            wind_kw=min(wind_kw, db.Const.wind_max_total_kw);
            for ib=1:numel(bess_mult)
                peak_vie=max(pv_kw,wind_kw); if peak_vie==0, peak_vie=base_opt; end
                bess_kw=min(round(peak_vie*bess_mult(ib)), db.Const.bess_max_kw);
                for ih=1:numel(bess_h)
                    bess_kwh=min(bess_kw*bess_h(ih), db.Const.bess_max_kwh);
                    ann_sav_o=0;
                    for si=1:4
                        P_PV_o=P_PV_s{si}*pv_kw;
                        P_Wind_o=P_Wind_s{si}*wind_kw;
                        [~,~,sav_s,~,~]=fast_dispatch(P_load_s{si},P_PV_o,P_Wind_o,...
                            bess_kw,bess_kwh,dg1_opt_val,dg2_opt_val,fp_o,et_o);
                        surplus_s=max(P_PV_o+P_Wind_o-P_load_s{si},0);
                        nm_rev_s=sum(surplus_s)*et_o*0.80;
                        ann_sav_o=ann_sav_o+(sav_s+nm_rev_s)*days_s(si);
                    end
                    % [FIX-C] Многоуровневый CAPEX
                    incr_capex_o = tiered_capex(pv_kw,   db.Const.capex_pv_tiers) ...
                                 + tiered_capex(wind_kw,  db.Const.capex_wind_tiers) ...
                                 + tiered_capex(bess_kwh, db.Const.capex_bess_tiers);
                    if ann_sav_o>0 && incr_capex_o>0
                        payback_o=incr_capex_o/ann_sav_o;
                        if payback_o<best_payback
                            best_payback=payback_o;
                            best_pv=pv_kw; best_wind=wind_kw;
                            best_bess_p=bess_kw; best_bess_e=bess_kwh;
                        end
                    end
                end
            end
        end
    end
    md=getappdata(fig,'modelData');
    md.pv_opt=best_pv; md.wind_opt=best_wind;
    md.bess_p_opt=best_bess_p; md.bess_e_opt=best_bess_e;
    md.dg1_opt=dg1_opt_val; md.dg2_opt=dg2_opt_val;
    md.dg_n_units=max(n1,n2);
    setappdata(fig,'modelData',md);
    pvLabel.Text=format_power(md.pv_opt);
    windLabel.Text=format_power(md.wind_opt);
    bessPLabel.Text=format_power(md.bess_p_opt);
    bessELabel.Text=format_energy(md.bess_e_opt);
    dg1Label.Text=format_power(md.dg1_opt);
    dg2Label.Text=format_power(md.dg2_opt);
    optBtn.Text='🔧 Авто-оптимизация ВИЭ';
    optBtn.BackgroundColor=[0.45 0.15 0.55];
end

% =========================================================================
% ОПТИМАЛЬНОЕ РАСПРЕДЕЛЕНИЕ НАГРУЗКИ МЕЖДУ ДГУ
% =========================================================================
function [P1_opt, P2_opt] = select_optimal_dispatch(Net, P1max, P2max, hours1, hours2, use1, use2)
    P1_opt=0; P2_opt=0;
    if Net<=0, return; end
    MIN_LOAD=0.30; SOFT_LO=0.50; SOFT_HI=0.85; FUEL_EPS=0.02;
    minP1=MIN_LOAD*P1max; minP2=MIN_LOAD*P2max;
    slo1=SOFT_LO*P1max; slo2=SOFT_LO*P2max;
    shi1=SOFT_HI*P1max; shi2=SOFT_HI*P2max;
    score_hours=@(p1,p2) p1*hours1+p2*hours2;

    function [sol,best_def] = try_dispatch(lo1,hi1,lo2,hi2)
        best_def=Inf; bf=Inf; bh=Inf; sol=[0,0];
        step=max(1,round(min(max(P1max,1),max(P2max,1))/80));
        if use1 && P1max>0
            P1c=min(max(Net,lo1),hi1);
            if P1c>=lo1-1e-6
                def=max(0,Net-P1c); fu=piecewise_fuel_consumption(P1c,P1max);
                hs=score_hours(P1c,0);
                [best_def,bf,bh,sol]=upd(def,fu,hs,[P1c,0],best_def,bf,bh,sol,FUEL_EPS);
            end
        end
        if use2 && P2max>0
            P2c=min(max(Net,lo2),hi2);
            if P2c>=lo2-1e-6
                def=max(0,Net-P2c); fu=piecewise_fuel_consumption(P2c,P2max);
                hs=score_hours(0,P2c);
                [best_def,bf,bh,sol]=upd(def,fu,hs,[0,P2c],best_def,bf,bh,sol,FUEL_EPS);
            end
        end
        if use1 && use2 && P1max>0 && P2max>0
            P1_lo_b=max(lo1,Net-hi2); P1_hi_b=min(hi1,Net-lo2);
            if P1_lo_b<=P1_hi_b+1e-6
                P1_list=unique([P1_lo_b:step:P1_hi_b, P1_hi_b]);
                for kk=1:length(P1_list)
                    P1c=P1_list(kk); P2c=Net-P1c;
                    if P2c>=lo2-1e-6 && P2c<=hi2+1e-6 && P1c>=lo1-1e-6
                        fu=piecewise_fuel_consumption(P1c,P1max)+piecewise_fuel_consumption(P2c,P2max);
                        hs=score_hours(P1c,P2c);
                        [best_def,bf,bh,sol]=upd(0,fu,hs,[P1c,P2c],best_def,bf,bh,sol,FUEL_EPS);
                    end
                end
            end
        end
    end

    [result,cur_def]=try_dispatch(slo1,shi1,slo2,shi2);
    if cur_def>1e-6
        [sol2,def2]=try_dispatch(minP1,shi1,minP2,shi2);
        if def2<cur_def-1e-6, result=sol2; cur_def=def2; end
    end
    if cur_def>1e-6
        [sol3,def3]=try_dispatch(minP1,P1max,minP2,P2max);
        if def3<cur_def-1e-6, result=sol3; end
    end
    P1_opt=result(1); P2_opt=result(2);
end

function [bd,bf,bh,bs]=upd(def,fuel,hs,sol,bd,bf,bh,bs,eps_rel)
    eps_f=max(eps_rel*max(fuel,1e-9),1e-9);
    if def<bd-1e-6
        bd=def; bf=fuel; bh=hs; bs=sol;
    elseif abs(def-bd)<=1e-6
        if fuel<bf-eps_f, bf=fuel; bh=hs; bs=sol;
        elseif abs(fuel-bf)<=eps_f && hs<bh-1e-3, bh=hs; bs=sol; end
    end
end

function [vie,lcoe,saving,def_pct,fuel_liters] = fast_dispatch(P_load,P_PV,P_Wind,bp,be,dg1c,dg2c,fp,et)
    eta_c=db.Const.eta_charge; eta_d=db.Const.eta_discharge;
    bess_self=db.Const.bess_self_disch;
    bess_reserve=0.20; DG_MIN=db.Const.dg_min_load_pct; DG_STAGE=db.Const.dg_stage_on_pct;
    curr_e=0.60*be;
    P_DG1=zeros(1,24); P_DG2=zeros(1,24);
    P_BESS=zeros(1,24); P_def=zeros(1,24);
    for h=1:24
        if be>0, curr_e=curr_e*(1-bess_self); end
        soc_h=curr_e/max(be,1)*100;
        Net=P_load(h)-P_PV(h)-P_Wind(h);
        if Net>0
            if soc_h>bess_reserve*100 && be>0
                max_disch=min(bp,(curr_e-bess_reserve*be)*eta_d);
                take=min(Net,max_disch);
                if take>0, P_BESS(h)=take; curr_e=curr_e-take/eta_d; Net=Net-take; end
            end
            if Net>0
                if dg1c>0 && dg2c>0
                    if Net>dg1c*DG_STAGE
                        total_c=dg1c+dg2c;
                        p1=min(dg1c,Net*dg1c/total_c); p2=min(dg2c,Net*dg2c/total_c);
                        leftover=Net-p1-p2;
                        p1_new=min(dg1c,p1+max(leftover,0));
                        leftover2=p1+max(leftover,0)-p1_new;
                        p2=min(dg2c,p2+max(leftover2,0));
                        P_DG1(h)=p1_new; P_DG2(h)=p2;
                    else
                        P1c=min(max(Net,DG_MIN*dg1c),dg1c);
                        P2c=min(max(Net,DG_MIN*dg2c),dg2c);
                        if piecewise_fuel_consumption(P1c,dg1c)<=piecewise_fuel_consumption(P2c,dg2c)
                            P_DG1(h)=min(Net,dg1c);
                        else
                            P_DG2(h)=min(Net,dg2c);
                        end
                    end
                    Net=max(0,Net-P_DG1(h)-P_DG2(h));
                elseif dg1c>0
                    P_DG1(h)=min(Net,dg1c); Net=Net-P_DG1(h);
                elseif dg2c>0
                    P_DG2(h)=min(Net,dg2c); Net=Net-P_DG2(h);
                end
            end
            if Net>0, P_def(h)=Net; end
        else
            surplus=-Net;
            if be>0
                max_charge=min(bp,(be-curr_e)/eta_c);
                charge=min(surplus,max_charge);
                if charge>0
                    curr_e=min(be,curr_e+charge*eta_c);
                    P_BESS(h)=-charge;
                end
            end
        end
    end
    E_load=sum(P_load); E_vie=sum(P_PV+P_Wind);
    E_served=E_load-sum(P_def);
    vie=min(100*E_vie/max(E_load,0.1),100);
    def_pct=100*sum(P_def)/max(E_load,0.1);
    fuel_liters=0;
    for h=1:24
        fuel_liters=fuel_liters+piecewise_fuel_consumption(P_DG1(h),dg1c)...
            +piecewise_fuel_consumption(P_DG2(h),dg2c);
    end
    % [FIX-DRY] OPEX из db.Const
    cost_total=fuel_liters*fp + sum(P_PV)*db.Const.opex_pv_per_kwh ...
        + sum(P_Wind)*db.Const.opex_wind_per_kwh + be*db.Const.opex_bess_per_kwh_yr/365;
    cost_grid=E_load*et;
    saving=cost_grid-cost_total;
    lcoe=cost_total/max(E_served,0.1);
end

% =========================================================================
% 8. ОСНОВНАЯ ФУНКЦИЯ СИМУЛЯЦИИ
% =========================================================================
function run_logic()
    t=0:23;
    base=max(loadField.Value,1);
    seed=max(1,sum(double(cityDrop.Value))+day(datePick.Value)+month(datePick.Value)*31+mod(year(datePick.Value),100)*365);
    rng(uint32(round(seed)));
    season=get_season(datePick.Value);
    day_type=get_day_type(datePick.Value);
    base_profile=db.Profiles.(season).(day_type);
    load_type_mult=db.LoadType.(loadTypeDrop.Value);
    base_profile=base_profile.*load_type_mult;
    prof_max=max(base_profile);
    if prof_max>0, base_profile=base_profile/prof_max; else, base_profile=ones(1,24); end
    holiday_k=holiday_multiplier(cityDrop.Value,datePick.Value);
    load_noise=1+ar1_noise(24,0.06,0.7);
    % P_load_det — всегда детерминированный реальный профиль нагрузки
    P_load_det=max(base*base_profile.*load_noise*holiday_k, 0);

    % [FIX-LSTM] LSTM используется ТОЛЬКО для управления BESS (bess_reserve),
    % а НЕ как замена реальной нагрузки.
    % Источник: Kong W. et al. IEEE Trans. Smart Grid 2019 — LSTM для EMS
    P_load=P_load_det;       % всегда реальный профиль
    P_load_lstm=P_load_det;  % прогноз для управления BESS
    lstm_used=false;
    try
        N_hist=length(allDaysData.P_load);
        if N_hist>=48
            P_prev_48h=allDaysData.P_load(end-47:end);
        elseif N_hist>0
            P_prev_48h=repmat(P_load_det,1,ceil(48/24));
            P_prev_48h=P_prev_48h(1:48);
        else
            P_prev_48h=repmat(P_load_det,1,2);
        end
        is_weekend=strcmp(get_day_type(datePick.Value),'Weekend');
        P_lstm=microgrid_lstm_predict(P_prev_48h,datePick.Value,is_weekend);
        if ~isempty(P_lstm) && length(P_lstm)==24 && all(isfinite(P_lstm))
            P_load_lstm=P_lstm;  % только для упреждающего управления BESS
            lstm_used=true;
        end
    catch
    end
    % P_load = P_load_det (реальный). P_load_lstm используется ниже только
    % для расчёта future_deficit при управлении BESS.

    md=getappdata(fig,'modelData');
    p_pv_val=md.pv_opt; p_wind_val=md.wind_opt;
    p_bess_val=md.bess_p_opt; e_bess=md.bess_e_opt;
    dg1_cap=md.dg1_opt; dg2_cap=md.dg2_opt;
    sources=md.sources;

    bess_soh=getappdata(fig,'bess_soh');
    if isempty(bess_soh), bess_soh=1.0; end
    e_bess_nom=e_bess; e_bess=e_bess*bess_soh;

    eta_c=db.Const.eta_charge; eta_d=db.Const.eta_discharge;
    bess_self=db.Const.bess_self_disch;
    ckey_reg=regDrop.Value;

    if sources.ses
        P_PV=calc_pv(t,p_pv_val,cityDrop.Value,datePick.Value,ckey_reg);
    else
        P_PV=zeros(1,24);
    end
    if sources.ves
        P_Wind=calc_wind(p_wind_val,season,t,ckey_reg);
    else
        P_Wind=zeros(1,24);
    end

    P_DG1=zeros(1,24); P_DG2=zeros(1,24);
    P_BESS=zeros(1,24); P_deficit=zeros(1,24); P_surplus=zeros(1,24);
    DG_MIN_LOAD_PCT=db.Const.dg_min_load_pct;

    last_soc=getappdata(fig,'last_soc_value');
    if isempty(last_soc), last_soc=80; end
    init_soc=last_soc;
    SOC=zeros(1,24); SOC(1)=init_soc;
    curr_e=max(0,min(e_bess,SOC(1)/100*e_bess));

    dg_hours=getappdata(fig,'dg_hours');
    hours1=dg_hours(1); hours2=dg_hours(2);

    for h=1:24
        if h>1 && e_bess>0
            curr_e=max(curr_e*(1-bess_self),0);
        end
        SOC(h)=max(0,min(100,(curr_e/max(e_bess,1))*100));

        % [FIX-LSTM] Управление BESS основано на прогнозе (P_load_lstm)
        h_remain=h:24;
        future_net=max(P_load_lstm(h_remain)-(P_PV(h_remain)+P_Wind(h_remain)),0);
        future_deficit=sum(future_net);
        future_peak=0;
        if length(h_remain)>1, future_peak=max(future_net(2:end)); end
        peak_guard=0;
        if e_bess>0 && future_peak>1.3*max(P_load(h),1), peak_guard=0.10; end
        if e_bess>0
            bess_reserve=0.20+0.15*(future_deficit/(e_bess+eps))+peak_guard;
            bess_reserve=min(bess_reserve,0.40);
        else
            bess_reserve=0.20;
        end

        % Реальная нагрузка P_load(h) — всегда P_load_det
        Net=P_load(h)-(P_PV(h)+P_Wind(h));

        if Net>0
            if SOC(h)>bess_reserve*100 && sources.bess && e_bess>0
                max_disch=min(p_bess_val,(curr_e-bess_reserve*e_bess)*eta_d);
                take=min(Net,max_disch);
                if take>0
                    P_BESS(h)=take; curr_e=curr_e-take/eta_d; Net=Net-take;
                end
            end
            need_dg1=sources.dgu1; need_dg2=sources.dgu2;
            if (need_dg1||need_dg2) && Net>0
                [P1_opt,P2_opt]=select_optimal_dispatch(Net,dg1_cap,dg2_cap,hours1,hours2,need_dg1,need_dg2);
                P_DG1(h)=P1_opt; P_DG2(h)=P2_opt; Net=Net-P1_opt-P2_opt;
                if P1_opt>0, hours1=hours1+1; end
                if P2_opt>0, hours2=hours2+1; end
            end
            % [FIX-2] Дозаряд BESS от ДГУ при минимальной нагрузке (IEC 62619:2022)
            if sources.bess && e_bess>0
                if need_dg1 && P_DG1(h)>0 && P_DG1(h)<dg1_cap*DG_MIN_LOAD_PCT && dg1_cap>0
                    excess=dg1_cap*DG_MIN_LOAD_PCT-P_DG1(h);
                    P_DG1(h)=dg1_cap*DG_MIN_LOAD_PCT;
                    max_charge=min(p_bess_val,(e_bess-curr_e)/eta_c);
                    actual_charge=min(excess,max_charge);
                    if actual_charge>0
                        P_BESS(h)=max(-p_bess_val,P_BESS(h)-actual_charge);
                        curr_e=min(e_bess,curr_e+actual_charge*eta_c);
                    end
                end
                if need_dg2 && P_DG2(h)>0 && P_DG2(h)<dg2_cap*DG_MIN_LOAD_PCT && dg2_cap>0
                    excess=dg2_cap*DG_MIN_LOAD_PCT-P_DG2(h);
                    P_DG2(h)=dg2_cap*DG_MIN_LOAD_PCT;
                    max_charge=min(p_bess_val,(e_bess-curr_e)/eta_c);
                    actual_charge=min(excess,max_charge);
                    if actual_charge>0
                        P_BESS(h)=max(-p_bess_val,P_BESS(h)-actual_charge);
                        curr_e=min(e_bess,curr_e+actual_charge*eta_c);
                    end
                end
            end
            Net=P_load(h)-(P_PV(h)+P_Wind(h)+max(P_BESS(h),0)+P_DG1(h)+P_DG2(h));
            if Net>0, P_deficit(h)=Net;
            elseif Net<-1e-6, P_surplus(h)=P_surplus(h)+(-Net); end
        else
            surplus=-Net;
            if sources.bess && e_bess>0
                max_charge_power=min(p_bess_val,(e_bess-curr_e)/eta_c);
                charge=min(surplus,max_charge_power);
                if charge>0
                    curr_e=min(e_bess,curr_e+charge*eta_c);
                    P_BESS(h)=-charge; surplus=surplus-charge;
                end
            end
            P_surplus(h)=max(surplus,0);
        end
    end

    if batch_mode, setappdata(fig,'dg_hours',[hours1,hours2]); end
    soc_final=max(0,min(100,(curr_e/max(e_bess,1))*100));
    setappdata(fig,'last_soc_value',soc_final);

    E_load=sum(P_load); E_deficit=sum(P_deficit);
    reg_key=regDrop.Value;
    if isfield(db.Econ,reg_key)
        re=db.Econ.(reg_key);
        fuel_price=get_city_fuel(cityDrop.Value,reg_key);
        elec_tariff=re.tariff; co2_grid=re.co2;
    else
        fuel_price=get_city_fuel(cityDrop.Value,reg_key);
        elec_tariff=8.5; co2_grid=0.40;
    end

    % [FIX-DRY] Все OPEX-константы из db.Const
    co2_kg_per_l     = db.Const.co2_kg_per_l_diesel;
    pv_opex          = db.Const.opex_pv_per_kwh;
    wind_opex        = db.Const.opex_wind_per_kwh;
    bess_opex_yr     = db.Const.opex_bess_per_kwh_yr;
    voll_rub_per_kwh = db.Const.voll_rub_per_kwh;

    E_PV=sum(P_PV); E_Wind=sum(P_Wind);
    E_surplus=sum(P_surplus); E_curtailed=E_surplus;
    grid_connected=false;

    fuel_liters=0;
    for hh=1:24
        fuel_liters=fuel_liters+piecewise_fuel_consumption(P_DG1(hh),dg1_cap)...
            +piecewise_fuel_consumption(P_DG2(hh),dg2_cap);
    end
    cost_fuel=fuel_liters*fuel_price;
    E_PV_useful=max(E_PV-E_curtailed*E_PV/(max(E_PV+E_Wind,1)),0);
    E_Wind_useful=max(E_Wind-E_curtailed*E_Wind/(max(E_PV+E_Wind,1)),0);
    cost_pv_op  =E_PV_useful*pv_opex;
    cost_wind_op=E_Wind_useful*wind_opex;
    if sources.bess && e_bess>0
        cost_bess_op=e_bess*bess_opex_yr/365;
    else
        cost_bess_op=0;
    end
    cost_total=cost_fuel+cost_pv_op+cost_wind_op+cost_bess_op;
    cost_grid_eq=(E_load-E_deficit)*elec_tariff;
    cost_saving=cost_grid_eq-cost_total;
    co2_DG_kg=fuel_liters*co2_kg_per_l;

    E_covered=max(E_load-E_deficit,1);
    % [FIX-E] LCOE включает стоимость недопоставки (VOLL)
    % Источник: ФАС РФ Методические рекомендации 2023; ENTSO-E VOLL Study 2022
    cost_with_voll = cost_total + E_deficit * voll_rub_per_kwh;
    LCOE=cost_with_voll/max(E_load,1);   % LCOE с учётом VOLL на всю нагрузку
    LCOE_no_voll=cost_total/E_covered;   % классический LCOE (для сравнения)
    vie_share=min((E_PV+E_Wind)/max(E_load,1)*100,100.0);

    curtail_loss=E_curtailed*LCOE_no_voll;
    if grid_connected
        net_metering_rate=elec_tariff*0.80;
        revenue_net_meter=E_surplus*net_metering_rate;
    else
        net_metering_rate=0; revenue_net_meter=0;
    end
    cost_saving_nm=cost_saving+revenue_net_meter;

    % [FIX-C] CAPEX с многоуровневым расчётом по масштабу
    capex_pv   = tiered_capex(p_pv_val,   db.Const.capex_pv_tiers);
    capex_wind = tiered_capex(p_wind_val,  db.Const.capex_wind_tiers);
    capex_bess = tiered_capex(e_bess_nom,  db.Const.capex_bess_tiers);
    capex_dg   = tiered_capex(dg1_cap,    db.Const.capex_dg_tiers) ...
               + tiered_capex(dg2_cap,    db.Const.capex_dg_tiers);
    capex_total=capex_pv+capex_wind+capex_bess+capex_dg;
    annual_saving_nm=cost_saving_nm*365;

    % [FIX-D] Финансовые параметры из db.Const (обновлено ЦБ РФ апрель 2026)
    discount_rate  = db.Const.discount_rate;   % 25% = 21% (ЦБ) + 4% риск
    project_life   = db.Const.project_life;
    fuel_inflation = db.Const.fuel_inflation;
    om_inflation   = db.Const.om_inflation;

    incremental_capex=capex_pv+capex_wind+capex_bess;

    % [FIX-7] bess_repl_year = 8 (IEC 62619:2022 ≈ 3000 цикл/365 = 8.2 г.)
    bess_replacement_year=db.Const.bess_repl_year;
    bess_replacement_capex=capex_bess*db.Const.bess_repl_cost_k;

    % [FIX-8] Остаточная стоимость с линейной амортизацией (IRENA 2023)
    pv_life=30; wind_life=25; dg_life=15; bess_life=bess_replacement_year;
    salvage_pv  =capex_pv  *max(0,(pv_life  -project_life)/pv_life);
    salvage_wind=capex_wind*max(0,(wind_life -project_life)/wind_life);
    salvage_dg  =capex_dg  *max(0,(dg_life   -project_life)/dg_life);
    bess_age_at_end=project_life-bess_replacement_year;
    salvage_bess=capex_bess*db.Const.bess_repl_cost_k*max(0,(bess_life-bess_age_at_end)/bess_life);
    salvage_value=salvage_pv+salvage_wind+salvage_dg+salvage_bess;

    diesel_fuel_liters_d=0;
    dg_cap_d=max(dg1_cap+dg2_cap,1);
    for hh_d=1:24
        p_d=min(P_load(hh_d),dg_cap_d);
        if p_d>0
            p1d=p_d*dg1_cap/dg_cap_d; p2d=p_d*dg2_cap/dg_cap_d;
            diesel_fuel_liters_d=diesel_fuel_liters_d...
                +piecewise_fuel_consumption(p1d,dg1_cap)...
                +piecewise_fuel_consumption(p2d,dg2_cap);
        end
    end
    diesel_opex_daily=diesel_fuel_liters_d*fuel_price;
    diesel_mgmt_daily=diesel_opex_daily*db.Const.mgmt_diesel_pct;
    diesel_total_daily=diesel_opex_daily+diesel_mgmt_daily;

    hybrid_mgmt_daily=cost_total*db.Const.mgmt_hybrid_pct;
    hybrid_total_daily=cost_total+hybrid_mgmt_daily;

    annual_diesel_base=diesel_total_daily*365;
    fuel_share_hybrid=cost_fuel/max(cost_total,1e-9);
    annual_hybrid_base=hybrid_total_daily*365;

    % [FIX-1] NPV: раздельная индексация топлива/O&M
    % Источник: IRENA "Renewable Energy Finance" (2020) — раздельная индексация
    npv=-incremental_capex;
    lcoe_lcc_cost=incremental_capex; lcoe_lcc_energy=0;
    for yr=1:project_life
        yr_diesel=annual_diesel_base*(1+fuel_inflation)^yr;
        yr_hybrid_fuel=annual_hybrid_base*fuel_share_hybrid*(1+fuel_inflation)^yr;
        yr_hybrid_om  =annual_hybrid_base*(1-fuel_share_hybrid)*(1+om_inflation)^yr;
        yr_hybrid=yr_hybrid_fuel+yr_hybrid_om;
        yr_sav=yr_diesel-yr_hybrid;
        npv=npv+yr_sav/(1+discount_rate)^yr;
        lcoe_lcc_cost=lcoe_lcc_cost+yr_hybrid/(1+discount_rate)^yr;
        E_yr=(E_load-E_deficit)*365;
        lcoe_lcc_energy=lcoe_lcc_energy+E_yr/(1+discount_rate)^yr;
        if yr==bess_replacement_year
            npv=npv-bess_replacement_capex/(1+discount_rate)^yr;
            lcoe_lcc_cost=lcoe_lcc_cost+bess_replacement_capex/(1+discount_rate)^yr;
        end
    end
    npv=npv+salvage_value/(1+discount_rate)^project_life;
    LCOE_LCC=iff(lcoe_lcc_energy>0, lcoe_lcc_cost/lcoe_lcc_energy, 0);

    annual_savings_incr=max((diesel_total_daily-hybrid_total_daily)*365,0);
    if annual_savings_incr>0, payback_simple=incremental_capex/annual_savings_incr;
    else, payback_simple=Inf; end
    if payback_simple<=7, payback_effect='✅ ВЫГОДНО';
    elseif payback_simple<=15, payback_effect='⚠️ ПРИЕМЛЕМО';
    else, payback_effect='❌ НЕВЫГОДНО'; end
    if isinf(payback_simple)||payback_simple>50
        payback_str=['∞ (убыток) ' payback_effect];
    elseif npv>0
        payback_str=sprintf('%.1f лет (NPV=+%.1f млн) %s',payback_simple,npv/1e6,payback_effect);
    else
        payback_str=sprintf('%.1f лет (NPV=−%.1f млн) %s',payback_simple,abs(npv)/1e6,payback_effect);
    end
    payback_years=payback_simple;

    % Диагностика
    H_DG_const=1.2;
    P_DG1_avg=mean(P_DG1); P_DG2_avg=mean(P_DG2);
    S_total=max(dg1_cap+dg2_cap+p_pv_val+p_wind_val,1);
    H_sys=H_DG_const*(P_DG1_avg+P_DG2_avg)/S_total;
    diag_lines={};
    off_list={};
    if ~sources.ses,  off_list{end+1}='СЭС'; end
    if ~sources.ves,  off_list{end+1}='ВЭС'; end
    if ~sources.dgu1, off_list{end+1}='ДГУ1'; end
    if ~sources.dgu2, off_list{end+1}='ДГУ2'; end
    if ~sources.bess, off_list{end+1}='BESS'; end
    if isempty(off_list), diag_lines{end+1}='OK: Все источники активны.';
    else, diag_lines{end+1}=['ОТКЛ: ' strjoin(off_list,', ')]; end
    total_deficit=sum(P_deficit);
    if total_deficit>0
        diag_lines{end+1}='─────────────────';
        diag_lines{end+1}=sprintf('⚠ ДЕФИЦИТ: %d ч.',sum(P_deficit>0));
        diag_lines{end+1}=sprintf('Макс: %.1f кВт',max(P_deficit));
        diag_lines{end+1}=sprintf('VOLL: %.0f руб/сут',E_deficit*voll_rub_per_kwh);
    else
        if ~isempty(off_list)
            diag_lines{end+1}='─────────────────';
            diag_lines{end+1}='✓ Дефицита нет.';
        else
            diag_lines{end+1}='✓ Система в норме.';
        end
    end
    diag_lines{end+1}='─────────────────';
    diag_lines{end+1}=sprintf('Топливо: %d руб/л',fuel_price);
    if E_curtailed>0.1
        diag_lines{end+1}=sprintf('Сброс ВИЭ: %.0f кВт·ч/сут',E_curtailed);
        diag_lines{end+1}=sprintf('Потери: -%.0f руб/сут',curtail_loss);
    end
    if cost_saving>=0, diag_lines{end+1}=sprintf('Экономия: +%.0f руб/сут',cost_saving);
    else, diag_lines{end+1}=sprintf('УБЫТОК: -%.0f руб/сут',abs(cost_saving)); end
    if lstm_used
        diag_lines{end+1}='─────────────────';
        diag_lines{end+1}='📡 LSTM: упрежд. BESS';
    end
    diag_lines{end+1}='─────────────────';
    diagArea.Value=diag_lines;

    deficit_key=round(total_deficit);
    if total_deficit>0 && deficit_key~=lastShownDeficit
        lastShownDeficit=deficit_key;
        uialert(fig,sprintf('Дефицит в %d ч!',sum(P_deficit>0)),'Предупреждение','Icon','warning');
    elseif total_deficit==0, lastShownDeficit=0; end

    % Сохранение результатов
    lastResults.t=t; lastResults.P_load=P_load;
    lastResults.P_PV=P_PV; lastResults.P_Wind=P_Wind;
    lastResults.P_DG1=P_DG1; lastResults.P_DG2=P_DG2;
    lastResults.P_BESS=P_BESS; lastResults.SOC=SOC;
    lastResults.P_deficit=P_deficit; lastResults.P_surplus=P_surplus;
    lastResults.city=cityDrop.Value;
    lastResults.p_pv_val=p_pv_val; lastResults.p_wind_val=p_wind_val;
    lastResults.p_bess_val=p_bess_val;
    lastResults.dg1_cap=dg1_cap; lastResults.dg2_cap=dg2_cap;
    lastResults.e_bess=e_bess_nom; lastResults.e_bess_eff=e_bess;
    lastResults.dateVal=datePick.Value; lastResults.base=base;
    lastResults.load_type=loadTypeDrop.Value;
    lastResults.season=season; lastResults.day_type=day_type;
    lastResults.econ.fuel_price=fuel_price;
    lastResults.econ.elec_tariff=elec_tariff;
    lastResults.econ.co2_grid=co2_grid;
    lastResults.econ.fuel_liters=fuel_liters;
    lastResults.econ.cost_fuel=cost_fuel;
    lastResults.econ.cost_pv_op=cost_pv_op;
    lastResults.econ.cost_wind_op=cost_wind_op;
    lastResults.econ.cost_bess_op=cost_bess_op;
    lastResults.econ.cost_total=cost_total;
    lastResults.econ.cost_grid_eq=cost_grid_eq;
    lastResults.econ.cost_saving=cost_saving;
    lastResults.econ.LCOE=LCOE;
    lastResults.econ.LCOE_no_voll=LCOE_no_voll;
    lastResults.econ.co2_DG_kg=co2_DG_kg;
    lastResults.econ.vie_share=vie_share;
    lastResults.econ.E_surplus=E_surplus;
    lastResults.econ.E_curtailed=E_curtailed;
    lastResults.econ.curtail_loss=curtail_loss;
    lastResults.econ.grid_connected=grid_connected;
    lastResults.econ.net_metering_rate=net_metering_rate;
    lastResults.econ.revenue_net_meter=revenue_net_meter;
    lastResults.econ.cost_saving_nm=cost_saving_nm;
    lastResults.econ.annual_saving_nm=annual_saving_nm;
    lastResults.econ.capex_pv=capex_pv;
    lastResults.econ.capex_wind=capex_wind;
    lastResults.econ.capex_bess=capex_bess;
    lastResults.econ.capex_dg=capex_dg;
    lastResults.econ.capex_total=capex_total;
    lastResults.econ.payback_years=payback_years;
    lastResults.econ.payback_str=payback_str;
    lastResults.econ.npv=npv;
    lastResults.econ.payback_effect=payback_effect;
    lastResults.econ.bess_soh=bess_soh;
    lastResults.econ.e_bess_eff=e_bess;
    lastResults.econ.LCOE_LCC=LCOE_LCC;
    lastResults.econ.salvage_value=salvage_value;
    lastResults.econ.om_inflation=om_inflation;
    lastResults.econ.H_sys=H_sys;
    lastResults.econ.voll_cost=E_deficit*voll_rub_per_kwh;
    lastResults.lstm_used=lstm_used;

    current_date=datePick.Value;
    if isempty(allDaysData.dates)
        allDaysData.P_load=P_load; allDaysData.P_PV=P_PV;
        allDaysData.P_Wind=P_Wind; allDaysData.P_DG1=P_DG1;
        allDaysData.P_DG2=P_DG2; allDaysData.P_BESS=P_BESS;
        allDaysData.SOC=SOC; allDaysData.P_deficit=P_deficit;
        allDaysData.P_surplus=P_surplus; allDaysData.dates={current_date};
    else
        date_idx=find(cellfun(@(d) isequal(d,current_date),allDaysData.dates),1);
        if ~isempty(date_idx)
            si=(date_idx-1)*24+1; ei=date_idx*24;
            allDaysData.P_load(si:ei)=P_load; allDaysData.P_PV(si:ei)=P_PV;
            allDaysData.P_Wind(si:ei)=P_Wind; allDaysData.P_DG1(si:ei)=P_DG1;
            allDaysData.P_DG2(si:ei)=P_DG2; allDaysData.P_BESS(si:ei)=P_BESS;
            allDaysData.SOC(si:ei)=SOC; allDaysData.P_deficit(si:ei)=P_deficit;
            allDaysData.P_surplus(si:ei)=P_surplus;
        else
            first_date=allDaysData.dates{1};
            if current_date<first_date
                allDaysData.P_load=[P_load,allDaysData.P_load];
                allDaysData.P_PV=[P_PV,allDaysData.P_PV];
                allDaysData.P_Wind=[P_Wind,allDaysData.P_Wind];
                allDaysData.P_DG1=[P_DG1,allDaysData.P_DG1];
                allDaysData.P_DG2=[P_DG2,allDaysData.P_DG2];
                allDaysData.P_BESS=[P_BESS,allDaysData.P_BESS];
                allDaysData.SOC=[SOC,allDaysData.SOC];
                allDaysData.P_deficit=[P_deficit,allDaysData.P_deficit];
                allDaysData.P_surplus=[P_surplus,allDaysData.P_surplus];
                allDaysData.dates=[{current_date},allDaysData.dates];
            else
                allDaysData.P_load=[allDaysData.P_load,P_load];
                allDaysData.P_PV=[allDaysData.P_PV,P_PV];
                allDaysData.P_Wind=[allDaysData.P_Wind,P_Wind];
                allDaysData.P_DG1=[allDaysData.P_DG1,P_DG1];
                allDaysData.P_DG2=[allDaysData.P_DG2,P_DG2];
                allDaysData.P_BESS=[allDaysData.P_BESS,P_BESS];
                allDaysData.SOC=[allDaysData.SOC,SOC];
                allDaysData.P_deficit=[allDaysData.P_deficit,P_deficit];
                allDaysData.P_surplus=[allDaysData.P_surplus,P_surplus];
                allDaysData.dates=[allDaysData.dates,{current_date}];
            end
        end
    end
    if ~batch_mode, update_graphs_with_new_data(); end
end

% =========================================================================
% ПЕРЕСЧЁТ ЭКОНОМИКИ (все константы из db.Const — принцип DRY)
% =========================================================================
function recalc_economics()
    if isempty(lastResults), return; end
    r=lastResults;
    E_PV=sum(r.P_PV); E_Wind=sum(r.P_Wind);
    E_load=sum(r.P_load); E_deficit=sum(r.P_deficit);
    fuel_liters=0;
    for hh=1:24
        fuel_liters=fuel_liters+piecewise_fuel_consumption(r.P_DG1(hh),r.dg1_cap)...
            +piecewise_fuel_consumption(r.P_DG2(hh),r.dg2_cap);
    end
    cost_fuel=fuel_liters*r.econ.fuel_price;
    % [FIX-DRY] Все OPEX из db.Const, не хардкоды
    cost_pv_op  =E_PV  *db.Const.opex_pv_per_kwh;
    cost_wind_op=E_Wind*db.Const.opex_wind_per_kwh;
    cost_bess_op=r.e_bess*db.Const.opex_bess_per_kwh_yr/365;
    cost_total  =cost_fuel+cost_pv_op+cost_wind_op+cost_bess_op;
    E_covered   =max(E_load-E_deficit,1);
    % [FIX-E] LCOE с VOLL
    voll_rub=db.Const.voll_rub_per_kwh;
    LCOE=(cost_total+E_deficit*voll_rub)/max(E_load,1);
    LCOE_no_voll=cost_total/E_covered;
    cost_grid_eq=(E_load-E_deficit)*r.econ.elec_tariff;
    cost_saving=cost_grid_eq-cost_total;
    % [FIX-DRY] CO2 из db.Const
    co2_DG_kg=fuel_liters*db.Const.co2_kg_per_l_diesel;
    vie_share=min((E_PV+E_Wind)/max(E_load,1)*100,100);

    r.econ.fuel_liters=fuel_liters; r.econ.cost_fuel=cost_fuel;
    r.econ.cost_pv_op=cost_pv_op; r.econ.cost_wind_op=cost_wind_op;
    r.econ.cost_bess_op=cost_bess_op;
    r.econ.cost_total=cost_total;
    r.econ.LCOE=LCOE; r.econ.LCOE_no_voll=LCOE_no_voll;
    r.econ.cost_saving=cost_saving;
    r.econ.co2_DG_kg=co2_DG_kg; r.econ.vie_share=vie_share;
    r.econ.cost_saving_nm=cost_saving+r.econ.revenue_net_meter;
    r.econ.annual_saving_nm=r.econ.cost_saving_nm*365;
    r.econ.voll_cost=E_deficit*voll_rub;

    % [FIX-DRY] [FIX-D] Финансовые параметры из db.Const
    discount_rate_rc  =db.Const.discount_rate;
    project_life_rc   =db.Const.project_life;
    fuel_inflation_rc =db.Const.fuel_inflation;
    om_inflation_rc   =db.Const.om_inflation;

    incremental_capex_rc=r.econ.capex_pv+r.econ.capex_wind+r.econ.capex_bess;
    bess_repl_year_rc  =db.Const.bess_repl_year;
    bess_repl_capex_rc =r.econ.capex_bess*db.Const.bess_repl_cost_k;

    [~,~,diesel_total_annual_rc,~,~,~,~]=calculate_diesel_mode(r);
    hybrid_mgmt_daily_rc =cost_total*db.Const.mgmt_hybrid_pct;
    hybrid_total_daily_rc=cost_total+hybrid_mgmt_daily_rc;
    hybrid_total_annual_rc=hybrid_total_daily_rc*365;
    annual_savings_rc=max(diesel_total_annual_rc-hybrid_total_annual_rc,0);

    % [FIX-1] Раздельная индексация топлива/O&M
    fuel_share_rc=cost_fuel/max(cost_total,1e-9);
    npv_rc=-incremental_capex_rc;
    for yr_rc=1:project_life_rc
        yr_sav_fuel=annual_savings_rc*fuel_share_rc*(1+fuel_inflation_rc)^yr_rc;
        yr_sav_om  =annual_savings_rc*(1-fuel_share_rc)*(1+om_inflation_rc)^yr_rc;
        yr_sav_rc  =yr_sav_fuel+yr_sav_om;
        npv_rc=npv_rc+yr_sav_rc/(1+discount_rate_rc)^yr_rc;
        if yr_rc==bess_repl_year_rc
            npv_rc=npv_rc-bess_repl_capex_rc/(1+discount_rate_rc)^yr_rc;
        end
    end
    if annual_savings_rc>0, r.econ.payback_years=incremental_capex_rc/annual_savings_rc;
    else, r.econ.payback_years=Inf; end
    if r.econ.payback_years<=7, payback_effect_rc='✅ ВЫГОДНО';
    elseif r.econ.payback_years<=15, payback_effect_rc='⚠️ ПРИЕМЛЕМО';
    else, payback_effect_rc='❌ НЕВЫГОДНО'; end
    if isinf(r.econ.payback_years)||r.econ.payback_years>50
        r.econ.payback_str=['∞ (убыток) ' payback_effect_rc];
    elseif npv_rc>0
        r.econ.payback_str=sprintf('%.1f лет (NPV=+%.1f млн) %s',r.econ.payback_years,npv_rc/1e6,payback_effect_rc);
    else
        r.econ.payback_str=sprintf('%.1f лет (NPV=−%.1f млн) %s',r.econ.payback_years,abs(npv_rc)/1e6,payback_effect_rc);
    end
    r.econ.npv=npv_rc; r.econ.payback_effect=payback_effect_rc;
    lastResults=r;
end

% =========================================================================
% ОБНОВЛЕНИЕ ГРАФИКОВ
% =========================================================================
function update_graphs_with_new_data()
    if isempty(lastResults)||isempty(allDaysData.dates), return; end
    r=lastResults;
    N_all=length(allDaysData.P_load); t_all=0:(N_all-1); view_win=24;
    cur_dn_v=datenum(datePick.Value); cur_didx=0;
    for dvi=1:length(allDaysData.dates)
        if abs(datenum(allDaysData.dates{dvi})-cur_dn_v)<0.5, cur_didx=dvi; break; end
    end
    if cur_didx>0
        x_start=max((cur_didx-1)*24-0.5,-0.5);
        x_end=min(x_start+view_win,N_all-0.5);
    else
        x_end=N_all-0.5; x_start=max(x_end-view_win,-0.5);
    end
    dragState.startXLim=[x_start,x_end];
    idx_visible=(t_all>=x_start)&(t_all<=x_end);

    cla(ax2); hold(ax2,'on');
    area(ax2,t_all,allDaysData.SOC,'FaceColor',[0.2 0.7 0.2],'FaceAlpha',0.7);
    plot(ax2,[t_all(1) t_all(end)],[20 20],'r--','LineWidth',1.5);
    plot(ax2,[t_all(1) t_all(end)],[80 80],'b--','LineWidth',1.0);
    hold(ax2,'off');
    title(ax2,'SOC батареи (%)'); ylim(ax2,[0 110]);
    legend(ax2,{'SOC','Мин 20%','Рек.80%'},'Location','best','FontSize',7);
    xlim(ax2,[x_start x_end]); apply_cyclic_xaxis(ax2,N_all);

    P_BESS_disch_all=max(allDaysData.P_BESS,0);
    P_BESS_charg_all=max(-allDaysData.P_BESS,0);
    P_curtail_all=zeros(size(allDaysData.P_load));
    if isfield(allDaysData,'P_surplus')&&~isempty(allDaysData.P_surplus)
        P_curtail_all=allDaysData.P_surplus;
    end

    visible_load=allDaysData.P_load(idx_visible);
    P_stk_all=allDaysData.P_PV+allDaysData.P_Wind+allDaysData.P_DG1+allDaysData.P_DG2+P_BESS_disch_all;
    visible_stk=P_stk_all(idx_visible);
    y_max_data=max([visible_load,visible_stk,P_curtail_all(idx_visible),allDaysData.P_deficit(idx_visible)]);
    if y_max_data<1, y_max_data=1; end
    y_top=y_max_data*1.18;
    mag=10^floor(log10(y_top)); y_nice=ceil(y_top/mag)*mag;

    cla(ax3); hold(ax3,'on');
    h_bars=bar(ax3,t_all,[allDaysData.P_PV;allDaysData.P_Wind;allDaysData.P_DG1;allDaysData.P_DG2;P_BESS_disch_all]','stacked');
    colors=[0.95 0.75 0.1;0.3 0.65 0.95;0.85 0.3 0.1;0.95 0.5 0.1;0.4 0.4 0.95];
    for c=1:5, h_bars(c).FaceColor=colors(c,:); end
    h_deficit=[]; h_curtail=[]; h_bess_chg=[];
    if any(allDaysData.P_deficit>0)
        h_deficit=bar(ax3,t_all,allDaysData.P_deficit,'FaceColor',[0.9 0.15 0.15],'EdgeColor',[0.7 0 0],'FaceAlpha',0.75);
    end
    if any(P_curtail_all>0.1)
        h_curtail=bar(ax3,t_all,P_curtail_all,'FaceColor',[0.6 0.6 0.6],'EdgeColor',[0.4 0.4 0.4],'FaceAlpha',0.45,'BarWidth',0.4);
    end
    if any(P_BESS_charg_all>0)
        h_bess_chg=plot(ax3,t_all,P_BESS_charg_all,'c--','LineWidth',1.5);
    end
    h_load=plot(ax3,t_all,allDaysData.P_load,'k-','LineWidth',2.5);
    ylim(ax3,[0 y_nice]); xlim(ax3,[x_start x_end]);
    apply_cyclic_xaxis(ax3,N_all); grid(ax3,'on');

    dg1_avg_pct=iff(r.dg1_cap>0,mean(r.P_DG1/r.dg1_cap)*100,0);
    dg2_avg_pct=iff(r.dg2_cap>0,mean(r.P_DG2/r.dg2_cap)*100,0);
    cur_md=getappdata(fig,'modelData');
    if ~isempty(cur_md), cur_src=cur_md.sources;
    else, cur_src=struct('ses',true,'ves',true,'dgu1',true,'dgu2',true,'bess',true); end
    dg1_str=iff(cur_src.dgu1,sprintf('ДГУ1:%s(%.0f%%)',format_power(r.dg1_cap),dg1_avg_pct),'ДГУ1:ОТКЛ');
    dg2_str=iff(cur_src.dgu2,sprintf('ДГУ2:%s(%.0f%%)',format_power(r.dg2_cap),dg2_avg_pct),'ДГУ2:ОТКЛ');
    title_line1=sprintf('%s [%s] | СЭС:%s ВЭС:%s %s %s BESS:%s/%s',...
        r.city,r.load_type,format_power(r.p_pv_val),format_power(r.p_wind_val),...
        dg1_str,dg2_str,format_power(r.p_bess_val),format_energy(r.e_bess));
    econ_label=iff(r.econ.cost_saving>=0,...
        sprintf('Экономия:+%.0f руб',r.econ.cost_saving),...
        sprintf('УБЫТОК:-%.0f руб',abs(r.econ.cost_saving)));
    curtail_str=iff(isfield(r.econ,'E_curtailed')&&r.econ.E_curtailed>0.1,...
        sprintf(' | Сброс:%.0fкВт·ч',r.econ.E_curtailed),'');
    title_line2=sprintf('ВИЭ:%.1f%% LCOE:%.2f руб CO₂:%.0f кг %s%s',...
        r.econ.vie_share,r.econ.LCOE,r.econ.co2_DG_kg,econ_label,curtail_str);
    title(ax3,{title_line1;title_line2},'FontSize',9);

    legend_handles=[]; legend_labels={};
    if any(allDaysData.P_PV>0.1),            legend_handles(end+1)=h_bars(1); legend_labels{end+1}='СЭС'; end
    if any(allDaysData.P_Wind>0.1),          legend_handles(end+1)=h_bars(2); legend_labels{end+1}='ВЭС'; end
    if cur_src.dgu1&&any(allDaysData.P_DG1>0.1), legend_handles(end+1)=h_bars(3); legend_labels{end+1}='ДГУ1'; end
    if cur_src.dgu2&&any(allDaysData.P_DG2>0.1), legend_handles(end+1)=h_bars(4); legend_labels{end+1}='ДГУ2'; end
    if any(P_BESS_disch_all>0.1),            legend_handles(end+1)=h_bars(5); legend_labels{end+1}='BESS↑'; end
    if ~isempty(h_deficit),                  legend_handles(end+1)=h_deficit;  legend_labels{end+1}='Дефицит'; end
    if ~isempty(h_curtail),                  legend_handles(end+1)=h_curtail;  legend_labels{end+1}='Сброс ВИЭ'; end
    if ~isempty(h_bess_chg),                 legend_handles(end+1)=h_bess_chg; legend_labels{end+1}='BESS↓заряд'; end
    legend_handles(end+1)=h_load; legend_labels{end+1}='Нагрузка';
    legend(ax3,legend_handles,legend_labels,'Location','northwest','FontSize',9);
    hold(ax3,'off');

    P_surplus_all=zeros(size(allDaysData.P_load));
    if isfield(allDaysData,'P_surplus')&&~isempty(allDaysData.P_surplus)
        P_surplus_all=allDaysData.P_surplus;
    end
    P_net_all=max(allDaysData.P_PV+allDaysData.P_Wind+allDaysData.P_DG1+...
        allDaysData.P_DG2+allDaysData.P_BESS-P_surplus_all, 0);
    cla(ax4); hold(ax4,'on');
    plot(ax4,t_all,P_net_all,'b-','LineWidth',2);
    plot(ax4,t_all,allDaysData.P_load,'k--','LineWidth',1.5);
    if any(allDaysData.P_deficit>0)
        area(ax4,t_all,allDaysData.P_deficit,'FaceColor',[1 0.5 0.5],'FaceAlpha',0.5);
    end
    title(ax4,'Генерация vs Нагрузка (с учётом BESS)'); grid(ax4,'on');
    legend(ax4,{'Нетто-генерация','Нагрузка','Дефицит'},'Location','best','FontSize',8);
    xlim(ax4,[x_start x_end]);
    visible_net=P_net_all(idx_visible);
    visible_load_ax4=allDaysData.P_load(idx_visible);
    y_max_ax4=max([visible_net,visible_load_ax4]);
    if y_max_ax4<1, y_max_ax4=1; end
    ylim(ax4,[0,ceil(y_max_ax4*1.05)]);
    apply_cyclic_xaxis(ax4,N_all); hold(ax4,'off');

    cla(ax1);
    plot(ax1,t_all,allDaysData.P_load,'-ro','MarkerFaceColor','r','LineWidth',1.5);
    title(ax1,'Нагрузка (кВт)'); grid(ax1,'on');
    xlim(ax1,[x_start x_end]); ylim(ax1,[0,max(visible_load)*1.05+1]);
    apply_cyclic_xaxis(ax1,N_all);

    cla(ax5);
    energy_pv=sum(r.P_PV); energy_wind=sum(r.P_Wind);
    energy_dgu1=sum(r.P_DG1); energy_dgu2=sum(r.P_DG2);
    energy_bess_out=sum(max(r.P_BESS,0));
    shares=[energy_pv,energy_wind,energy_dgu1,energy_dgu2,energy_bess_out];
    total_energy=sum(shares);
    base_labels={'СЭС','ВЭС','ДГУ1','ДГУ2','BESS'};
    labels_with_percent=cell(1,length(shares));
    for i=1:length(shares)
        if total_energy>0
            labels_with_percent{i}=sprintf('%s (%.1f%%)',base_labels{i},(shares(i)/total_energy)*100);
        else
            labels_with_percent{i}=base_labels{i};
        end
    end
    pie_colors=[0.95 0.75 0.1;0.3 0.65 0.95;0.85 0.3 0.1;0.95 0.5 0.1;0.4 0.4 0.95];
    mask=shares>0.1;
    if any(mask)
        h_pie=pie(ax5,shares(mask),labels_with_percent(mask));
        idx=find(mask); patch_count=0;
        for i=1:length(h_pie)
            if isa(h_pie(i),'matlab.graphics.primitive.Patch')
                patch_count=patch_count+1;
                if patch_count<=length(idx), h_pie(i).FaceColor=pie_colors(idx(patch_count),:); end
            end
            if isa(h_pie(i),'matlab.graphics.primitive.Text')
                h_pie(i).FontSize=11; h_pie(i).FontWeight='bold';
            end
        end
    else
        text(ax5,0,0,'Нет данных','HorizontalAlignment','center','FontSize',12);
    end
    drawnow;
end

% =========================================================================
% 9. ОТЧЁТЫ
% =========================================================================
function open_report()
    if isempty(lastResults)
        uialert(fig,'Сначала запустите моделирование!','Нет данных','Icon','warning'); return;
    end
    r=lastResults;
    rfig=uifigure('Name',sprintf('Сравнительный отчёт — %s %s',r.city,string(r.dateVal,'dd.MM.yyyy')),...
        'Position',[50 50 1400 850]); rfig.Color=[1 1 1];
    uilabel(rfig,'Position',[20 800 1200 30],...
        'Text',sprintf('СРАВНИТЕЛЬНЫЙ ЭКОНОМИЧЕСКИЙ ОТЧЁТ — %s — %s',r.city,string(r.dateVal,'dd.MM.yyyy')),...
        'FontSize',16,'FontWeight','bold','FontColor',[0.1 0.3 0.6]);
    uibutton(rfig,'Position',[1180 795 200 35],'Text','🖨️ Печать отчёта',...
        'BackgroundColor',[0.3 0.5 0.8],'FontColor','w','FontWeight','bold','FontSize',12,...
        'ButtonPushedFcn',@(b,e) print_report(rfig));
    uilabel(rfig,'Position',[20 770 1360 20],...
        'Text',sprintf('Нагрузка: %s | Дата: %s | Тип: %s | Ставка: %.0f%% | VOLL: %.0f руб/кВт·ч',...
            format_power(r.base),string(r.dateVal,'dd.MM.yyyy'),r.load_type,...
            db.Const.discount_rate*100, db.Const.voll_rub_per_kwh),...
        'FontSize',10,'FontColor',[0.3 0.3 0.3]);
    uipanel(rfig,'Position',[20 760 1360 2],'BackgroundColor',[0.5 0.5 0.5]);

    capex_hybrid=r.econ.capex_pv+r.econ.capex_wind+r.econ.capex_bess+r.econ.capex_dg;
    capex_diesel=r.econ.capex_dg;
    additional_capex=capex_hybrid-capex_diesel;
    capex_data={
        'СЭС',r.econ.capex_pv/1e6,0,r.econ.capex_pv/1e6;
        'ВЭС',r.econ.capex_wind/1e6,0,r.econ.capex_wind/1e6;
        'BESS',r.econ.capex_bess/1e6,0,r.econ.capex_bess/1e6;
        'ДГУ1+ДГУ2',r.econ.capex_dg/1e6,r.econ.capex_dg/1e6,0;
        'Итого CAPEX',capex_hybrid/1e6,capex_diesel/1e6,additional_capex/1e6};
    uitable(rfig,'Position',[20 640 1360 110],...
        'ColumnName',{'Компонент','Гибрид (млн руб)','Только ДГУ (млн руб)','Доп. CAPEX (млн руб)'},...
        'Data',capex_data,'FontSize',10);

    [~,~,diesel_total_cost,diesel_management,diesel_dg_pct,diesel_co2,diesel_fuel]=calculate_diesel_mode(r);
    [~,~,hybrid_total_cost,hybrid_management,hybrid_dg_pct,hybrid_co2,hybrid_fuel]=calculate_hybrid_mode(r);
    diesel_fuel_opex=diesel_fuel*r.econ.fuel_price;
    opex_data={
        'Топливо ДГУ',r.econ.cost_fuel*365/1e3,diesel_fuel_opex/1e3,(diesel_fuel_opex-r.econ.cost_fuel*365)/1e3;
        'Обслуживание СЭС',r.econ.cost_pv_op*365/1e3,0,r.econ.cost_pv_op*365/1e3;
        'Обслуживание ВЭС',r.econ.cost_wind_op*365/1e3,0,r.econ.cost_wind_op*365/1e3;
        'Обслуживание BESS',r.econ.cost_bess_op*365/1e3,0,r.econ.cost_bess_op*365/1e3;
        'Управл. расходы *',hybrid_management/1e3,diesel_management/1e3,(diesel_management-hybrid_management)/1e3;
        'VOLL (недопост.)',r.econ.voll_cost*365/1e3,0,0;
        'Итого OPEX',hybrid_total_cost/1e3,diesel_total_cost/1e3,(diesel_total_cost-hybrid_total_cost)/1e3};
    uitable(rfig,'Position',[20 490 1360 140],...
        'ColumnName',{'Статья','Гибрид (тыс руб/год)','Только ДГУ (тыс руб/год)','Экономия (тыс руб/год)'},...
        'Data',opex_data,'FontSize',10);
    uilabel(rfig,'Position',[20 480 1360 14],...
        'Text','* Управленческие расходы: гибрид 12% OPEX, только ДГУ 15% OPEX. Ставка: ЦБ РФ 21%+4% риск=25% (апрель 2026).',...
        'FontSize',8,'FontColor',[0.5 0.5 0.5]);

    if sum(r.P_load)>0, lcoe_diesel=diesel_total_cost/(sum(r.P_load)*365/1e3);
    else, lcoe_diesel=0; end
    npv_value=r.econ.npv/1e6;
    npv_str=iff(npv_value>=0,sprintf('+%.2f',npv_value),sprintf('−%.2f',abs(npv_value)));
    pb_eff=iff(isfield(r.econ,'payback_effect'),r.econ.payback_effect,'—');
    summary_data={
        'Нагрузка базовая',format_power(r.base),'—','—';
        'Средняя загрузка ДГУ (%)',sprintf('%.1f%%',hybrid_dg_pct),sprintf('%.1f%%',diesel_dg_pct),sprintf('−%.1f п.п.',diesel_dg_pct-hybrid_dg_pct);
        'Выбросы CO₂ (т/год)',sprintf('%.1f',hybrid_co2),sprintf('%.1f',diesel_co2),sprintf('−%.1f т/год (%.0f%%)',diesel_co2-hybrid_co2,iff(diesel_co2>0,(diesel_co2-hybrid_co2)/diesel_co2*100,0));
        'Расход топлива (тыс л/год)',sprintf('%.1f',hybrid_fuel/1000),sprintf('%.1f',diesel_fuel/1000),sprintf('−%.1f тыс.л/год',(diesel_fuel-hybrid_fuel)/1000);
        'Доля ВИЭ (%)',sprintf('%.1f%%',r.econ.vie_share),'0%',sprintf('+%.1f%%',r.econ.vie_share);
        'LCOE с VOLL (руб/кВт·ч)',sprintf('%.2f',r.econ.LCOE),sprintf('%.2f',lcoe_diesel),'—';
        'LCOE без VOLL (руб/кВт·ч)',sprintf('%.2f',r.econ.LCOE_no_voll),'—','—';
        'Срок окупаемости',r.econ.payback_str,'—',pb_eff;
        'NPV (млн руб)',npv_str,'—','—';
        'Ставка дисконтир.',sprintf('%.0f%%',db.Const.discount_rate*100),'—','ЦБ РФ 21%+4% риск'};
    uitable(rfig,'Position',[20 310 1360 165],...
        'ColumnName',{'Показатель','Гибрид','Только ДГУ','Изменение'},...
        'Data',summary_data,'FontSize',10);
    uilabel(rfig,'Position',[20 20 1360 18],...
        'Text','Microgrid AI v22.2 | IEC 62898-1:2017, ISO 8528-1:2018, IEC 61400-1:2019, IRENA 2023, BloombergNEF 2024, ФАС РФ VOLL 2023',...
        'FontSize',8,'FontColor',[0.6 0.6 0.6],'HorizontalAlignment','center');
end

function open_dgu_report()
    if isempty(lastResults)
        uialert(fig,'Сначала запустите моделирование!','Нет данных','Icon','warning'); return;
    end
    r=lastResults;
    dguFig=uifigure('Name',sprintf('Отчёт по ДГУ — %s %s',r.city,string(r.dateVal,'dd.MM.yyyy')),...
        'Position',[100 100 1400 850]); dguFig.Color=[1 1 1];
    uilabel(dguFig,'Position',[20 800 1200 30],...
        'Text',sprintf('ДЕТАЛЬНЫЙ ОТЧЁТ ПО ДГУ — %s — %s',r.city,string(r.dateVal,'dd.MM.yyyy')),...
        'FontSize',16,'FontWeight','bold','FontColor',[0.1 0.3 0.6]);
    uilabel(dguFig,'Position',[20 770 1360 20],...
        'Text',sprintf('Базовая нагрузка: %s | СЭС: %s | ВЭС: %s | BESS: %s / %s',...
            format_power(r.base),format_power(r.p_pv_val),format_power(r.p_wind_val),...
            format_power(r.p_bess_val),format_energy(r.e_bess)),...
        'FontSize',10,'FontColor',[0.3 0.3 0.3]);
    uipanel(dguFig,'Position',[20 760 1360 2],'BackgroundColor',[0.5 0.5 0.5]);

    ax_dgu=uiaxes(dguFig,'Position',[50 580 1300 160]);
    hours=0:23; dg1_pct=zeros(1,24); dg2_pct=zeros(1,24);
    for h=1:24
        if r.dg1_cap>0, dg1_pct(h)=(r.P_DG1(h)/r.dg1_cap)*100; end
        if r.dg2_cap>0, dg2_pct(h)=(r.P_DG2(h)/r.dg2_cap)*100; end
    end
    hold(ax_dgu,'on');
    bar(ax_dgu,hours,[dg1_pct;dg2_pct]','grouped');
    ylabel(ax_dgu,'Загрузка (%)'); xlabel(ax_dgu,'Час суток');
    title(ax_dgu,'Почасовая загрузка ДГУ (в % от номинала единичного агрегата)','FontSize',12,'FontWeight','bold');
    legend(ax_dgu,{sprintf('ДГУ1 (%s)',format_power(r.dg1_cap)),sprintf('ДГУ2 (%s)',format_power(r.dg2_cap))},'Location','northeast');
    grid(ax_dgu,'on'); ylim(ax_dgu,[0 110]); xlim(ax_dgu,[-0.5 23.5]);
    xticks(ax_dgu,0:2:23); hold(ax_dgu,'off');

    fuel_per_hour=zeros(1,24);
    for h=1:24
        fuel_per_hour(h)=piecewise_fuel_consumption(r.P_DG1(h),r.dg1_cap)...
            +piecewise_fuel_consumption(r.P_DG2(h),r.dg2_cap);
    end
    dg_cap_total=max(r.dg1_cap+r.dg2_cap,1);
    tableData=cell(24,7);
    for h=1:24
        tableData{h,1}=sprintf('%d:00',h-1);
        tableData{h,2}=sprintf('%.1f',r.P_DG1(h));
        tableData{h,3}=sprintf('%.1f%%',dg1_pct(h));
        tableData{h,4}=sprintf('%.1f',r.P_DG2(h));
        tableData{h,5}=sprintf('%.1f%%',dg2_pct(h));
        tableData{h,6}=sprintf('%.2f',fuel_per_hour(h));
        tableData{h,7}=sprintf('%.1f',(r.P_DG1(h)+r.P_DG2(h))/dg_cap_total*100);
    end
    uitable(dguFig,'Position',[20 220 1360 340],...
        'ColumnName',{'Час','ДГУ1 (кВт)','ДГУ1 (%)','ДГУ2 (кВт)','ДГУ2 (%)','Расход (л/ч)','Общ. загрузка (%)'},...
        'ColumnWidth',{60,90,80,90,80,90,100},'Data',tableData);

    total_fuel=sum(fuel_per_hour);
    total_energy_dg=sum(r.P_DG1+r.P_DG2);
    co2_kg=total_fuel*db.Const.co2_kg_per_l_diesel;
    dg_share=iff(sum(r.P_load)>0,sum(r.P_DG1+r.P_DG2)/sum(r.P_load)*100,0);

    stats_panel=uipanel(dguFig,'Position',[20 70 1360 140],'Title','Сводная статистика','FontSize',12,'FontWeight','bold');
    uilabel(stats_panel,'Position',[10 95 450 25],'Text',sprintf('ДГУ1: сред.%.1f%%, макс.%.1f%%, работа %dч',mean(dg1_pct),max(dg1_pct),sum(r.P_DG1>0)),'FontSize',11);
    uilabel(stats_panel,'Position',[10 65 450 25],'Text',sprintf('ДГУ2: сред.%.1f%%, макс.%.1f%%, работа %dч',mean(dg2_pct),max(dg2_pct),sum(r.P_DG2>0)),'FontSize',11);
    uilabel(stats_panel,'Position',[10 35 400 25],'Text',sprintf('Общая загрузка (сред.): %.1f%%',mean((r.P_DG1+r.P_DG2)/dg_cap_total*100)),'FontSize',11);
    uilabel(stats_panel,'Position',[480 95 430 25],'Text',sprintf('Выработка: %.0f кВт·ч/сут (%.0f тыс.кВт·ч/год)',total_energy_dg,total_energy_dg*365/1000),'FontSize',11);
    uilabel(stats_panel,'Position',[480 65 430 25],'Text',sprintf('Расход: %.1f л/сут (%.1f тыс.л/год)',total_fuel,total_fuel*365/1000),'FontSize',11);
    uilabel(stats_panel,'Position',[480 35 430 25],'Text',sprintf('Выбросы CO₂: %.1f кг/сут (%.1f т/год)',co2_kg,co2_kg*365/1000),'FontSize',11);
    uilabel(stats_panel,'Position',[940 95 380 25],'Text',sprintf('Доля ДГУ: %.1f%% | Топливо: %.0f руб/сут',dg_share,total_fuel*r.econ.fuel_price),'FontSize',11);
    uilabel(stats_panel,'Position',[940 65 380 25],'Text',sprintf('Ед. мощность ДГУ1: %s | ДГУ2: %s',format_power(r.dg1_cap),format_power(r.dg2_cap)),'FontSize',11);

    uibutton(dguFig,'Position',[1200 20 180 35],'Text','🖨️ Печать отчёта',...
        'BackgroundColor',[0.3 0.5 0.8],'FontColor','w','FontWeight','bold',...
        'ButtonPushedFcn',@(b,e) print_report(dguFig));
end

function open_dgu_chars()
    if isempty(lastResults)
        uialert(fig,'Сначала запустите моделирование!','Нет данных','Icon','warning'); return;
    end
    r=lastResults;
    P1_nom=r.dg1_cap; P2_nom=r.dg2_cap;
    if P1_nom==0&&P2_nom==0
        uialert(fig,'Оба ДГУ отключены!','Нет данных','Icon','warning'); return;
    end
    if P1_nom>0
        P_range1=0:1:P1_nom;
        B1=arrayfun(@(p) piecewise_fuel_consumption(p,P1_nom),P_range1);
        dBdP1=gradient(B1,P_range1); alpha1_deg=atand(dBdP1);
    else
        P_range1=0; B1=0; alpha1_deg=0;
    end
    if P2_nom>0
        P_range2=0:1:P2_nom;
        B2=arrayfun(@(p) piecewise_fuel_consumption(p,P2_nom),P_range2);
        dBdP2=gradient(B2,P_range2); alpha2_deg=atand(dBdP2);
    else
        P_range2=0; B2=0; alpha2_deg=0;
    end
    charFig=uifigure('Name','Характеристики ДГУ: B(P) и α(P)','Position',[200 200 1100 750]);
    charFig.Color=[1 1 1];
    axB=uiaxes(charFig,'Position',[60 380 980 300]);
    hold(axB,'on');
    if P1_nom>0
        plot(axB,P_range1,B1,'b-','LineWidth',2,'DisplayName',sprintf('ДГУ1 (%s)',format_power(P1_nom)));
        idx70=find(P_range1>=0.7*P1_nom,1,'first');
        if idx70>1&&idx70<length(P_range1)
            x_tan=0.7*P1_nom+[-30,30];
            y_tan=B1(idx70)+dBdP1(idx70)*(x_tan-0.7*P1_nom);
            plot(axB,x_tan,y_tan,'b--','LineWidth',1.2);
            text(0.7*P1_nom+5,B1(idx70)+5,sprintf('α=%.1f°',alpha1_deg(idx70)),'FontSize',9,'Color','b');
        end
    end
    if P2_nom>0
        plot(axB,P_range2,B2,'r-','LineWidth',2,'DisplayName',sprintf('ДГУ2 (%s)',format_power(P2_nom)));
        idx70_2=find(P_range2>=0.7*P2_nom,1,'first');
        if idx70_2>1&&idx70_2<length(P_range2)
            x_tan2=0.7*P2_nom+[-30,30];
            y_tan2=B2(idx70_2)+dBdP2(idx70_2)*(x_tan2-0.7*P2_nom);
            plot(axB,x_tan2,y_tan2,'r--','LineWidth',1.2);
            text(0.7*P2_nom+5,B2(idx70_2)+5,sprintf('α=%.1f°',alpha2_deg(idx70_2)),'FontSize',9,'Color','r');
        end
    end
    hold(axB,'off');
    xlabel(axB,'Мощность P, кВт'); ylabel(axB,'Расход B, л/ч');
    title(axB,'Расходная характеристика ДГУ B(P)'); grid(axB,'on'); legend(axB,'Location','northwest');

    axAlpha=uiaxes(charFig,'Position',[60 50 980 280]);
    hold(axAlpha,'on');
    if P1_nom>0, plot(axAlpha,P_range1,alpha1_deg,'b-','LineWidth',2,'DisplayName',sprintf('ДГУ1 (%s)',format_power(P1_nom))); end
    if P2_nom>0, plot(axAlpha,P_range2,alpha2_deg,'r-','LineWidth',2,'DisplayName',sprintf('ДГУ2 (%s)',format_power(P2_nom))); end
    hold(axAlpha,'off');
    xlabel(axAlpha,'Мощность P, кВт'); ylabel(axAlpha,'Угол α, градусы');
    title(axAlpha,'Угол α(P) = atan(dB/dP)'); grid(axAlpha,'on'); legend(axAlpha,'Location','northeast');
    uilabel(charFig,'Position',[60 20 980 25],...
        'Text','α = arctan(dB/dP). Касательные при 70% нагрузки. ISO 8528-1:2018 — аппроксимация кусочно-линейной характеристикой.',...
        'FontSize',10,'FontColor',[0.4 0.4 0.4],'HorizontalAlignment','center');
end

% =========================================================================
% ВСПОМОГАТЕЛЬНЫЕ ФУНКЦИИ
% =========================================================================
function B = piecewise_fuel_consumption(P, Pnom)
    % Кусочно-линейная аппроксимация характеристики ДГУ
    % Источник: ISO 8528-5:2013 — Performance and measurement
    % Пять отрезков: 0-30%, 30-50%, 50-75%, 75-85%, 85-100% нагрузки
    if P<=0||Pnom<=0, B=0; return; end
    B0  =0.04*0.28*Pnom;
    B_30=B0  +0.380*(0.30*Pnom);
    B_50=B_30+0.290*(0.20*Pnom);
    B_75=B_50+0.265*(0.25*Pnom);
    B_85=B_75+0.285*(0.10*Pnom);
    p=P/Pnom;
    if     p<=0.30, B=B0  +0.380*P;
    elseif p<=0.50, B=B_30+0.290*(P-0.30*Pnom);
    elseif p<=0.75, B=B_50+0.265*(P-0.50*Pnom);
    elseif p<=0.85, B=B_75+0.285*(P-0.75*Pnom);
    else,           B=B_85+0.360*(P-0.85*Pnom);
    end
end

function res = iff(condition, true_val, false_val)
    if condition, res=true_val; else, res=false_val; end
end

function [capex,opex,total_cost,management,dg_pct,co2,fuel] = calculate_diesel_mode(r)
    dg1_cap=r.dg1_cap; dg2_cap=r.dg2_cap;
    % [FIX-C] CAPEX ДГУ с многоуровневым расчётом
    capex=tiered_capex(dg1_cap,db.Const.capex_dg_tiers)+tiered_capex(dg2_cap,db.Const.capex_dg_tiers);
    P_dg_total=zeros(1,24);
    dg_cap_total=max(dg1_cap+dg2_cap,1);
    for h=1:24, P_dg_total(h)=min(r.P_load(h),dg_cap_total); end
    dg_pct=mean(P_dg_total/dg_cap_total)*100;
    fuel=0;
    for h=1:24
        p_h=P_dg_total(h);
        if p_h>0
            p1=p_h*dg1_cap/dg_cap_total; p2=p_h*dg2_cap/dg_cap_total;
            fuel=fuel+piecewise_fuel_consumption(p1,dg1_cap)+piecewise_fuel_consumption(p2,dg2_cap);
        end
    end
    fuel=fuel*365;
    co2=fuel*db.Const.co2_kg_per_l_diesel/1000;
    opex=fuel*r.econ.fuel_price;
    management=opex*db.Const.mgmt_diesel_pct;
    total_cost=opex+management;
end

function [capex,opex,total_cost,management,dg_pct,co2,fuel] = calculate_hybrid_mode(r)
    capex=r.econ.capex_total;
    bess_op=iff(isfield(r.econ,'cost_bess_op'),r.econ.cost_bess_op,0);
    opex=r.econ.cost_fuel*365+(r.econ.cost_pv_op+r.econ.cost_wind_op+bess_op)*365;
    management=opex*db.Const.mgmt_hybrid_pct;
    total_cost=opex+management;
    dg_total_cap=r.dg1_cap+r.dg2_cap;
    dg_pct=iff(dg_total_cap>0,(mean(r.P_DG1+r.P_DG2)/dg_total_cap)*100,0);
    co2=r.econ.co2_DG_kg*365/1000;
    fuel=r.econ.fuel_liters*365;
end

function print_report(target_fig)
    temp_pdf=[tempname '.pdf'];
    try
        drawnow;
        if exist('exportapp','file')==2
            exportapp(target_fig,temp_pdf);
        else
            error('exportapp недоступен');
        end
    catch
        try
            frame=getframe(target_fig);
            f_temp=figure('Visible','off','Position',[100 100 1000 800]);
            imshow(frame.cdata); print(f_temp,temp_pdf,'-dpdf','-r300'); close(f_temp);
        catch
            uialert(target_fig,'Не удалось экспортировать отчёт.','Ошибка','Icon','error'); return;
        end
    end
    try
        if ispc, winopen(temp_pdf);
        elseif ismac, system(['open "' temp_pdf '"']);
        else, system(['xdg-open "' temp_pdf '"']); end
        uialert(target_fig,sprintf('Отчёт сохранён:\n%s',temp_pdf),'Печать','Icon','success');
    catch
        uialert(target_fig,'Файл создан, но не удалось открыть автоматически.','Info','Icon','info');
    end
end

% =========================================================================
% 10. НАВИГАЦИЯ ПО ГРАФИКУ
% =========================================================================
function apply_cyclic_xaxis(ax, N)
    if N==0, return; end
    x_lo=ax.XLim(1); x_hi=ax.XLim(2);
    view_range=max(x_hi-x_lo,1);
    step=max(1,ceil(view_range/12));
    first_tick=ceil(x_lo/step)*step;
    hours_in_view=first_tick:step:floor(x_hi);
    hours_in_view=hours_in_view(hours_in_view>=0&hours_in_view<N);
    if isempty(hours_in_view), return; end
    labels=cell(1,length(hours_in_view));
    for kk=1:length(hours_in_view)
        labels{kk}=sprintf('%02d:00',mod(hours_in_view(kk),24));
    end
    ax.XTick=hours_in_view; ax.XTickLabel=labels;
    if ~isempty(allDaysData.dates)
        n_days=length(allDaysData.dates);
        first_day_idx=max(1,floor(hours_in_view(1)/24)+1);
        last_day_idx=min(n_days,floor(hours_in_view(end)/24)+1);
        if first_day_idx==last_day_idx
            date_str=datestr(allDaysData.dates{first_day_idx},'dd.mm.yyyy');
        else
            date_str=sprintf('%s … %s',...
                datestr(allDaysData.dates{first_day_idx},'dd.mm'),...
                datestr(allDaysData.dates{last_day_idx},'dd.mm'));
        end
        xlabel(ax,sprintf('Время суток (часы) — %s',date_str),'FontSize',9);
    else
        xlabel(ax,'Время суток (часы)','FontSize',9);
    end
end

function on_drag_start(~,~)
    cp=fig.CurrentObject; if isempty(cp), return; end
    parent_ax=cp;
    while ~isempty(parent_ax)&&~isa(parent_ax,'matlab.ui.control.UIAxes')
        parent_ax=parent_ax.Parent;
    end
    if isempty(parent_ax)||~(parent_ax==ax1||parent_ax==ax2||parent_ax==ax3||parent_ax==ax4), return; end
    pos=fig.CurrentPoint;
    dragState.active=true; dragState.startPixX=pos(1);
    dragState.startXLim=ax3.XLim; dragState.dir=0;
end

function on_drag_move(~,~)
    if ~dragState.active, return; end
    N=length(allDaysData.P_load); if N==0, return; end
    pos=fig.CurrentPoint; dpix=pos(1)-dragState.startPixX;
    dragState.dir=dpix;
    ax_w=ax3.Position(3); orig_r=dragState.startXLim(2)-dragState.startXLim(1);
    dx_dat=-dpix/max(ax_w,1)*orig_r;
    new_s=max(dragState.startXLim(1)+dx_dat,-0.5);
    new_e=new_s+orig_r;
    if new_e>N-0.5, new_e=N-0.5; new_s=max(new_e-orig_r,-0.5); end
    lim=[new_s,new_e];
    xlim(ax1,lim); xlim(ax2,lim); xlim(ax3,lim); xlim(ax4,lim);
end

function on_drag_end(~,~)
    if ~dragState.active, return; end
    dragState.active=false;
    N=length(allDaysData.P_load); if N==0, return; end
    at_left_edge=(ax3.XLim(1)<=0.5);
    at_right_edge=(ax3.XLim(2)>=N-0.6);
    swiped_left=(dragState.dir<-30); swiped_right=(dragState.dir>30);
    if swiped_left&&at_right_edge
        datePick.Value=datePick.Value+caldays(1);
        old_cb=datePick.ValueChangedFcn; datePick.ValueChangedFcn=[];
        run_logic(); datePick.ValueChangedFcn=old_cb;
    end
    if swiped_right&&at_left_edge
        datePick.Value=datePick.Value-caldays(1);
        old_cb=datePick.ValueChangedFcn; datePick.ValueChangedFcn=[];
        run_logic(); datePick.ValueChangedFcn=old_cb;
    end
    apply_cyclic_xaxis(ax1,length(allDaysData.P_load));
    apply_cyclic_xaxis(ax2,length(allDaysData.P_load));
    apply_cyclic_xaxis(ax3,length(allDaysData.P_load));
    apply_cyclic_xaxis(ax4,length(allDaysData.P_load));
end

% =========================================================================
% ЗАПУСК
% =========================================================================
run_logic();

fig.CloseRequestFcn=@(src,evt) on_fig_close(src);
function on_fig_close(src)
    if isappdata(fig,'last_soc_value'), rmappdata(fig,'last_soc_value'); end
    if isappdata(0,'MicrogridFigure'), rmappdata(0,'MicrogridFigure'); end
    delete(src);
end
end