Does anyone know how to modify this program for peak load shaving using charging and discharging battery

I am tryint to modify this code to shave the profile_load by energy storage battery (charging / discharging). The figure below this is my goal. Where is the orange curve is load_profile and the blue curve is charging/ discharging battery. The data in the figure below is different from the my load_profile.

Any help would be appreciated.

*********************************************************************************************************************************************************

load_profile=[250;227;206;220;197;203;210;191;353;464;518;544;588;565;562;606;576;494;409;322;301;288;268;228];

loss=0.01;

T=24;

M=1;

ESS_Cmax=90;

ESS_ch = optimvar('ESS_ch',T,'LowerBound',0,'UpperBound',M);

ESS_disch = optimvar('ESS_disch',T,'LowerBound',0,'UpperBound',M);

ESS_SOC = optimvar('ESS_SOC',T,'LowerBound',0,'UpperBound',ESS_Cmax);

Grid_E = optimvar('Grid_E',T,'LowerBound',0);

prob.Constraints.energyStorage = optimconstr(T);

prob.Constraints.energyStorage(1) = ESS_SOC(1) == 0;

prob.Constraints.energyStorage(2:T) = ESS_SOC(2:T) == ESS_SOC(1:T-1)*(1-loss) + ESS_ch(1:T-1) - ESS_disch(1:T-1);

prob.Constraints.EnergyBalance = Grid_E == load_profile - ESS_ch + ESS_disch;

%create binary variables

Z_X = optimvar('Z_X',T,'Type','integer','LowerBound',0,'UpperBound',1); % charging binary variable array

Z_Y = optimvar('Z_Y',T,'Type','integer','LowerBound',0,'UpperBound',1); % discharging binary variable array

% Consraints for binary variables

prob.Constraints.Ch_on = ESS_ch <= M*Z_X;

prob.Constraints.Dis_on = ESS_disch <= M*Z_Y;

prob.Constraints.notPositive = Z_X + Z_Y <= 1;

%**************************************************

W = optimvar('W');

Z = optimvar('Z');

prob = optimproblem;

prob.Objective = W - Z + 1/100*(sum(Z_X + Z_Y));

prob.Constraints.shaveW = W >= Grid_E;

prob.Constraints.shaveZ = Z <= Grid_E;

prob.ObjectiveSense = 'minimize';

options = optimoptions(prob.optimoptions,'Display','none')

[sol,fval,exitflag,output,lambda] = solve(prob,'Options',options)

% **************************************************

[Alan Weiss]

You set M = 1. Why? I suggest that you try 100. Or whatever your maximum value really is.

If you don't understand the code that you are using, feel free to ask what various specific lines mean. It appears to me that you are using https://www.mathworks.com/matlabcentral/answers/579006-objective-function-for-peak-shaving-and-valley-filling-of-wind-turbine-electrical-load-profile?s_tid=srchtitle .

Alan Weiss

MATLAB mathematical toolbox documentation