Error using .* : Conversion to optim.prob​lemdef.Opt​imizationE​xpression from optim.prob​lemdef.Opt​imizationI​nequality is not possible

Error using .*

Conversion to optim.problemdef.OptimizationExpression from optim.problemdef.OptimizationInequality is not possible.

Error in power_opti (line 78)

E_bat=cumsum(delta_t*P_bat.*(P_bat<=0)/eta_d_bat + delta_t*P_bat.*(P_bat>=0)/eta_c_bat)

I used optimproblem object to solve optimization problme but i end up with having this error.

%%%%%whole code

close all;
clear all;
clc;
load('wave_p.mat');
%%%%%%55555  non variable input
P_owc=Wave_KW;;%%%%% [10000x1]
P_load= 40+10*rand(length(Wave_KW),1)    ;%%%%% [10000x1]

wavepower_optimization=optimproblem;

delta_t=0.0224;
%%%% variables 
% HESS
% Total power output from the HESS will be the sum of power capacity of each ESS 

P_hess=P_owc-P_load;
%P_hess(t) = P_bat(t) + P_fw(t) + P_sc(t);
P_bat=zeros(length(Wave_KW),1);
P_fw=zeros(length(Wave_KW),1);
P_sc=zeros(length(Wave_KW),1);

%%%%%%%%%
P_bat_lower=0;
P_fw_lower=0;
P_sc_lower=0;

P_bat_upper=100;
P_fw_upper=100;
P_sc_upper=100;

r=0.06;
y_fw=15;
y_sc=25;
l_owc=10;
u_sc=1;
u_bat=1;
u_fw=1;
m_fw=0.9;
m_bat=0.7;
m_sc=0.9;


eta_d_bat=0.8;
eta_c_bat=0.8;

eta_d_fw=0.8;
eta_c_fw=0.8;

eta_d_bat=0.8;
eta_c_bat=0.8;

eta_d_sc=0.8;
eta_c_sc=0.8;


E_bat_min=0;
E_bat_max=1200;
E_fw_min=0;
E_fw_max=1200;

E_sc_min=0;
E_sc_max=1200;

%%%%%%%555
P_bat=optimvar('P_bat',length(Wave_KW),'LowerBound',P_bat_lower,'UpperBound',P_bat_upper);
P_fw=optimvar('P_fw',length(Wave_KW),'LowerBound',P_fw_lower,'UpperBound',P_fw_upper);
P_sc=optimvar('P_sc',length(Wave_KW),'LowerBound',P_sc_lower,'UpperBound',P_sc_upper);

%%%% p cosntrains
wavepower_optimization.Constraints.econs1=P_bat+P_fw+P_sc==P_hess
%%%%%% define


% Battery pack
%Discharging mode;
%E_bat(t+1)=E_bat(t)-(delta_t.*P_bat(t))/eta_d_bat;  t = 1, 2,….T;

E_bat=cumsum(delta_t*P_bat.*(P_bat<=0)/eta_d_bat + delta_t*P_bat.*(P_bat>=0)/eta_c_bat)
%E_bat=func_E_bat(delta_t,P_bat,eta_d_bat,eta_c_bat);

%Charging mode;
%E_bat(t+1)=E_bat(t)+ (delta_t.*P_bat(t).*eta_c_bat; t = 1, 2,….T;

%P_bat = E_bat/delta_t;  % battery power
SOC_bat= E_bat./E_bat_max;  % battery SOC

% Flywheel
%Discharging mode;
%E_fw(t+1)=E_fw(t)-(delta_t.*P_fw(t))/eta_d_fw;  t = 1, 2,….T;
%Charging mode;
%E_fw(t+1)=E_fw(t)+ (delta_t.*P_fw(t).eta_c_fw; t = 1, 2,….T;
E_fw=cumsum(delta_t*P_fw.*(P_fw<=0)/eta_d_fw + delta_t*P_fw.*(P_fw>=0)/eta_c_fw)

%E_fw=func_E_fw(delta_t,P_fw,eta_d_fw,eta_c_fw);

%P_fw = E_fw/delta_t;  % Flywheel power
SOC_fw= E_fw./E_fw_max;  % Flywheel SOC


% Supercapacitor bank
%Discharging mode;
%E_sc(t+1)=E_sc(t)-(delta_t.*P_sc(t))/eta_d_sc;  t = 1, 2,….T;
%Charging mode;
%E_sc(t+1)=E_sc(t)+ (delta_t.*P_sc(t).eta_c_sc; t = 1, 2,….T;
E_sc=cumsum(delta_t.*P_sc.*(P_sc<=0)./eta_d_sc + delta_t.*P_sc.*(P_sc>=0)./eta_c_sc)
%E_sc=func_E_sc(delta_t,P_sc,eta_d_sc,eta_c_sc)


%P_sc = E_sc/delta_t;  % SC power
SOC_sc= E_sc/E_sc_max;  % SC SOC

% Constraints
wavepower_optimization.Constraints.cons1=E_bat>=E_bat_min;
wavepower_optimization.Constraints.cons2=E_bat<=E_bat_max;

wavepower_optimization.Constraints.cons3=E_fw>=E_fw_min;
wavepower_optimization.Constraints.cons4=E_fw<=E_fw_max;


wavepower_optimization.Constraints.cons5=E_sc>=E_sc_min;
wavepower_optimization.Constraints.cons6=E_sc<=E_sc_max;
% Power balance
%P_owc(t)+ P_hess(t) = P_load(t)

% OWC constraints 
% % turbine generated power range
%omega_owc_min <= omega_owc(t) <= omega_owc_max % turbine safe operation speed range
%%%% why
% HESS constraints
%P_hess_min <= P_hess(t) <= P_hess_max
%E_hess_min <= E_hess(t) <= E_hess_max
%SOC_hess_min <= SOC_hess(t) <= SOC_hess_max
%%%%why

% Battery constraints
%P_bat_min <= P_bat(t) <= P_bat_max
%E_bat_min <= E_bat(t) <= E_bat_max
%SOC_bat_min <= SOC_bat(t) <= SOC_bat_max


% Flywheel constraints
%P_fw_min <= P_fw(t) <= P_fw_max
%E_fw_min <= E_fw(t) <=E_fw_max 
%SOC_fw_min <= SOC_fw(t) <= SOC_fw_max
%omega_fw_min <= omega_fw(t) <= omega_fw_max



% SC constraints
%P_sc_min <= P_sc(t) <= P_sc_max
%E_sc_min <= E_sc(t) <= E_sc_max
%SOC_sc_min <= SOC_sc(t) <= SOC_sc_max
% Loss of power supply probability (LPSP)
%0 < LPSP <= LPSP_desired

%P_supply(t).* delta_t = (P_owc(t).*delta_t + E_hess(t-1) - E_hess_min).*eta_inv
%LPS(t) = (P_load(t)-P_supply(t)).*delta_t when P_load(t) > P_supply(t)
%LPS(t) = 0 when P_load(t) < P_supply(t)
 %LPS  = func_LPS(delta_t,P_owc,P_load)
 
 LPS = (P_load-P_supply).*delta_t.*(P_load>P_supply)
%%%%%%%%%%
LPSP=sum(LPS./(P_load.*delta_t))*100 %? %cannot type the exact equation
%LPSP = func_LPSP(LPS,P_owc,P_load,delta_t)

wavepower_optimization.Constraints.cons2=LPSP<=LPSP_desired;

% NPC = Cap(cost) + O&M(cost)+ Rep(cost)
% Cap(Cost)=u_s*P_s*CRF where, u_s = Unit cost in $/kW, 

%Cap_Cost = func_cap_cost(u_s,P_sc,CRF)
% P_s = Power capacity in kW, s is the index for bat, fw and sc
% "Capital recovery factor" CRF = r*(1+r)^y/((1+r)^y-1) where r is the interest rate and y is the life span of each component 


% O&M(Cost)= m_s*u_s* P_s where m_s is a pre-defined percentage 
% Rep(Cost)= u_s*P_s*(l_owc-y)/y  where l_owc is life span of OWC plant and y is life span of each component 
   %battery


%NPC_bat = func_NPC_bat(P_bat,u_bat,y_bat,r,m_bat,l_owc);
CRF_bat = r*(1+r)^y_bat/((1+r)^y_bat-1);
NPC_bat = u_bat*P_bat*CRF_bat + m_bat*u_bat*P_bat + u_bat*P_bat*(l_owc-y_bat)/y_bat ; 

%NPC_sc= func_NPC_sc(P_sc,u_sc,y_sc,r,m_sc,l_owc);
CRF_fw = r*(1+r)^y_fw/((1+r)^y_fw-1);
NPC_fw = u_fw*P_fw*CRF_fw + m_fw*u_fw*P_fw+ u_fw*P_fw*(l_owc-y_fw)/y_fw  ;

%NPC_fw = func_NPC_fw(P_fw,u_fw,y_fw,r,m_fw,l_owc);
CRF_sc = r*(1+r)^y_sc/((1+r)^y_sc-1);
NPC_sc = u_sc*P_sc*CRF_sc + m_sc*u_sc*P_sc+ u_sc*P_sc*(l_owc-y_sc)/y_sc; 

%TCE=func_TCE(NPC_bat, NPC_sc,NPC_fw,P_load);
% Objective Function
%Min(TCE)

% Objective Function


TCE=(NPC_bat + NPC_sc + NPC_fw);

wavepower_optimization.Objective=TCE;
%%%%%%%%%%%%%%%%
x0.P_bat=randn(length(Wave_KW),1)*P_bat_max;
x0.P_fw=randn(length(Wave_KW),1)*P_fw_max;
x0.P_sc=randn(length(Wave_KW),1)*P_sc_max;

%x = particleswarm(problem)
[sol,mincost]=solve(wavepower_optimization,x0)

Answer 1

[Walter Roberson]

P_bat.*(P_bat<=0)

A relationship operator between two optimization expressions returns a constraint expression. Constraint expressions do not return binary values and cannot be used arithmetically such as by multiplying them by an optimization expression.

Unfortunately for you, min() and max() are not supported on optimization expressions either.

You are going to need to find a different approach.

[Walter Roberson]

(x-sqrt(x.^2))/2

is x for negative x, and 0 for positive x. Use + to get x for positive x and 0 for negative x.

These are supported operations for optimization expressions, but it might be the case that they are too complex to be supported when put together.

[kasun prabhath]

thank you. it works on this problem.

Answer 2

[Samira Mohammed Osman]

Error in optim.problemdef.OptimizationExpression/createBinary

[Samira Mohammed Osman]

clear all
close all 
%x_ij ∈ {0,1} ; Y_i ∈ {0,1} ∀ i,j
x11 = optimvar('x11','Type','integer','LowerBound',0, 'UpperBound',1);
x12 = optimvar('x21','Type','integer','LowerBound',0, 'UpperBound',1);
x13 = optimvar('x13','Type','integer','LowerBound',0, 'UpperBound',1);
x21 = optimvar('x21','Type','integer','LowerBound',0, 'UpperBound',1);
x22 = optimvar('x22','Type','integer','LowerBound',0, 'UpperBound',1);
x23 = optimvar('x23','Type','integer','LowerBound',0, 'UpperBound',1);
x31 = optimvar('x31','Type','integer','LowerBound',0, 'UpperBound',1);
x32 = optimvar('x32','Type','integer','LowerBound',0, 'UpperBound',1);
x33 = optimvar('x33','Type','integer','LowerBound',0, 'UpperBound',1);
Y1  = optimvar('Y1','Type','integer','LowerBound',0, 'UpperBound',1);
Y2  = optimvar('Y2','Type','integer','LowerBound',0, 'UpperBound',1);
Y3  = optimvar('Y3','Type','integer','LowerBound',0, 'UpperBound',1);
%Objective Function
prob = optimproblem('Objective',600*11*x11 + 500*12*x12+ 400*22*x13 + 600*8*x21+ 500*14*x22+ 400*10*x23 + 600*20*x31+ 500*10*x32+ 400*10*x33,'ObjectiveSense','min');

Error using optim.internal.problemdef.ExpressionTree/createBinary
OptimizationVariables appearing in the same OptimizationExpression must have distinct "Name" properties.
Make a new variable with a different "Name" property, or retrieve the original variable using the Variables property.

Error in optim.internal.problemdef.ExpressionForest/createBinary

Error in optim.problemdef.OptimizationExpression/createBinary

Error in +

%Constraints
prob.Constraints.c1 = x11 + x21 + x31 == 1;
prob.Constraints.c2 = x12 + x22 + x32 == 1;
prob.Constraints.c3 = x13 + x23 + x33 == 1;
prob.Constraints.c4 = Y1 + Y2 + Y3 == 2;
prob.Constraints.c5 = x11 <= Y1;
prob.Constraints.c6 = x12 <= Y1;
prob.Constraints.c7 = x13 <= Y1;
prob.Constraints.c8 = x21 <= Y2;
prob.Constraints.c9 = x22 <= Y2;
prob.Constraints.c10 = x23 <= Y2;
prob.Constraints.c11 = x31 <= Y3;
prob.Constraints.c12 = x32 <= Y3;
prob.Constraints.c13 = x33 <= Y3;
%convert
problem = prob2struct(prob);
[sol,fval,exitflag,output] = intlinprog(problem)

[Walter Roberson]

Notice your line

x12 = optimvar('x21','Type','integer','LowerBound',0, 'UpperBound',1);

You have x12 on the left side, but x21 on the right side. You also have

x21 = optimvar('x21','Type','integer','LowerBound',0, 'UpperBound',1);

which again uses x21 .