How to solve absolute problem in optimization

Hello, i want to solve a quadratic optimization problem (prob.Objective =sum((PgridV).^2);), in constrains i have a one variable where i need to find the absolute value of (Pb1dc). so i introduce a variable called 'K' . where -K<=Pb1dc<=K;

• FYI: i have used absoulte value in the following constrains (prob.Constraints.loadBalanceAC=Pb1==Pb1dc-(0.05*K);)

MATLAB solution showing exactely whati want . its working well.

However, when i want to chage the objective function with little modifcation (prob.Objective =sum((PgridV-M).^2);) where M is a reference signal (mean of Pload), when i run the simulation

absolute vaules (K) are not getting exactley what mean for. its showing random values ( K = absolute ((Pb1dc))

clc
clear all
Pload=[0;1;3;4;2;6;9;10;2;4]; % load 
N=10;
M=mean(Pload)+zeros(N,1); 
Einit1=0.5; % initial  energy 
E=zeros(N,1);
Emin1 = 0;       % mini energy 
Emax1 = 3;
dt=1;

prob = optimproblem;

PgridV = optimvar('PgridV',N,'LowerBound',0,'UpperBound',20); % grid power 
Pb1= optimvar('Pb1',N,'LowerBound',-1,'UpperBound',1); % ac power 
Pb1dc= optimvar('Pb1dc',N,'LowerBound',-1,'UpperBound',1); % dc power 
K=optimvar('K',N,'LowerBound',0);% absolute of (Pb1dc)
EbattV1 = optimvar('EbattV1',N,'LowerBound',Emin1,'UpperBound',Emax1); % energy 

prob.ObjectiveSense = 'minimize';
% prob.Objective =sum(K.^2);
%prob.Objective =sum((PgridV).^2);
prob.Objective =sum((PgridV-M).^2);

% 1 constrains
prob.Constraints.energyBalance = optimconstr(N);
prob.Constraints.energyBalance(1) = EbattV1(1) == Einit1-Pb1dc(1)*dt; %Its Ploss is constanat 
prob.Constraints.energyBalance(2:N-1) = EbattV1(2:N-1) == EbattV1(1:N-2)-Pb1dc(2:N-1)*dt;
prob.Constraints.energyBalance(N) = EbattV1(N) ==Einit1;

% K constrain for Pb1dc modulous 
prob.Constraints.kbalance1=optimconstr(N);
prob.Constraints.kbalance1(1:N)=-K(1:N)<=Pb1dc(1:N);
prob.Constraints.kbalance2=optimconstr(N);
prob.Constraints.kbalance2(1:N)=Pb1dc(1:N)<=K(1:N);

% load Balance
prob.Constraints.loadBalance=PgridV+Pb1==Pload;

% loss term
prob.Constraints.loadBalanceAC=Pb1==Pb1dc-(0.05*K);

options = optimoptions(prob.optimoptions,'Display','final');
% options = optimoptions(prob.optimoptions,'Algorithm','interior-point');
[values,fval,exitflag] = solve(prob,'Options',options)

Solving problem using lsqlin.

Minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the value of the optimality tolerance,
and constraints are satisfied to within the value of the constraint tolerance.

values = struct with fields:
    EbattV1: [10×1 double]
          K: [10×1 double]
        Pb1: [10×1 double]
      Pb1dc: [10×1 double]
     PgridV: [10×1 double]

fval = 57.4575

exitflag = 
    OptimalSolution

% Parse optmization results
if exitflag <= 0
    PgridV = zeros(N,1);
    Pb1 = zeros(N,1);
    Pb1dc = zeros(N,1);
    EbattV1 = zeros(N,1);
    K = zeros(N,1);

else
    PgridV = values.PgridV ;
    Pb1 =   values.Pb1;
    Pb1dc = values.Pb1dc
    EbattV1 = values.EbattV1;
    K = values.K    
end 

Pb1dc = 10×1

-0.6434 -0.6204 -0.6701 0.1224 -0.6886 1.0000 1.0000 1.0000 -0.8624 -0.0018

K = 10×1

7.1327 7.5914 6.5986 4.4490 6.2283 1.0000 1.0000 1.0000 2.7522 1.9640

.

[Torsten]

Can you show in your code above exactly what you mean ?

I ran the code with

prob.Objective =sum((PgridV-M).^2);

instead of

prob.Objective =sum((PgridV).^2);

(see above). What is not satisfactory with the results ?

[Jayachandra Malavatu]

[Jayachandra Malavatu]

The K values : which are absolute values of (Pb1dc)

[Torsten]

What's wrong with the K values listed above ?

They are not defined as the absolute values of Pb1dc in your constraints.

Your constraint is abs(Pb1dc) <= K.

Answer 1

[Jayachandra Malavatu]

[Jayachandra Malavatu]

Could you please help me how can write constrain if i want a variable K, which is abs(Pb1dc)

[Torsten]

Yes, it's easily seen from the output of your program (see above) that the constraints are satisfied:

abs(Pb1dc) <= K

[Torsten]

fun = @(x,y) x-abs(y);
prob.Constraints.kbalance1=optimconstr(N);
prob.Constraints.kbalance1(1:N)=fcn2optimexpr(fun,K,Pb1dc)==0;

[Torsten]

Didn't you read my answer ? By using "fcn2optimexpr" as done above. But note that this will make your problem nonlinear. Do you really need the absolute values of Pb1dc ? It will make your simulation much more complicated (although "abs" sounds so easy).

[Jayachandra Malavatu]

I have modified

fun = @(x,y) x-abs(y);

prob.Constraints.kbalance1=optimconstr(N);

prob.Constraints.kbalance1(1:N)=fcn2optimexpr(fun,K,Pb1dc)==0;

with

prob.Constraints.kbalance1=optimconstr(N);

prob.Constraints.kbalance1(1:N)=-K(1:N)<=Pb1dc(1:N);

prob.Constraints.kbalance2=optimconstr(N);

prob.Constraints.kbalance2(1:N)=Pb1dc(1:N)<=K(1:N);..

..

.. i got following error since i never used fmincon

Solving problem using fmincon.

Error using optim.problemdef.OptimizationProblem/solve

SOLVE requires a non-empty initial point structure to solve a nonlinear problem.

Error in modfun_effi_loss (line 48)

[values,~,exitflag] = solve(prob,'Options',options);

[Torsten]

Yes, add an initial guess for your solution variables in the call to "solve".

See the second example under

https://uk.mathworks.com/help/optim/ug/optim.problemdef.optimizationproblem.solve.html

on how to do this.

[Jayachandra Malavatu]

Yes, i need absolute value. because i tried initially with mixed integer , and my optimization problem become MINL, where i am unable to solve the issue with ga, (which i already tried), To overcome mixed integer i used this method.

[Torsten]

Without using fcn2optimexpr, you could add the constraints

K.^2 - Pb1dc.^2 == 0

with K constrained to be >=0.

[Jayachandra Malavatu]

Thanks . In my case i have a variables PgridV,Pb1,Pb1dc,K,EbattV1,

did i need to mention initial values for all above five variables with (N,1) like

x0.PgridV=3+zeros(N,1);

x0.Pb1=1+zeros(N,1);

x0.Pb1dc=1+zeros(N,1);

x0.K=1+zeros(N,1);

x0.EbattV1=1.5+zeros(N,1);

and finally at end

[values,~,exitflag] = solve(prob,x0,'Options',options);

is this something correct ?

[Torsten]

is this something correct ?

Yes. Is "lsqnonlin" chosen as solver ?

[Jayachandra Malavatu]

its fmincon

[Jayachandra Malavatu]

Is there any other way i can improve simulation time , since i need to consider N= 96 samples and optimizer taking a lot of time.

[Torsten]

clc
clear all
Pload=[0;1;3;4;2;6;9;10;2;4]; % load 
N=10;
M=mean(Pload)+zeros(N,1); 
Einit1=0.5; % initial  energy 
E=zeros(N,1);
Emin1 = 0;       % mini energy 
Emax1 = 3;
dt=1;

prob = optimproblem;

PgridV = optimvar('PgridV',N,'LowerBound',0,'UpperBound',20); % grid power 
Pb1= optimvar('Pb1',N,'LowerBound',-1,'UpperBound',1); % ac power 
Pb1dc= optimvar('Pb1dc',N,'LowerBound',-1,'UpperBound',1); % dc power 
K=optimvar('K',N,'LowerBound',0);% absolute of (Pb1dc)
EbattV1 = optimvar('EbattV1',N,'LowerBound',Emin1,'UpperBound',Emax1); % energy 

prob.ObjectiveSense = 'minimize';
% prob.Objective =sum(K.^2);
%prob.Objective =sum((PgridV).^2);
prob.Objective =sum((PgridV-M).^2);

% 1 constrains
prob.Constraints.energyBalance = optimconstr(N);
prob.Constraints.energyBalance(1) = EbattV1(1) == Einit1-Pb1dc(1)*dt; %Its Ploss is constanat 
prob.Constraints.energyBalance(2:N-1) = EbattV1(2:N-1) == EbattV1(1:N-2)-Pb1dc(2:N-1)*dt;
prob.Constraints.energyBalance(N) = EbattV1(N) ==Einit1;

% K constrain for Pb1dc modulous 
prob.Constraints.kbalance1=optimconstr(N);
prob.Constraints.kbalance1=K.^2-Pb1dc.^2==0;

% load Balance
prob.Constraints.loadBalance=PgridV+Pb1==Pload;

% loss term
prob.Constraints.loadBalanceAC=Pb1==Pb1dc-(0.05*K);

x0.PgridV=3+zeros(N,1);
x0.Pb1=1+zeros(N,1);
x0.Pb1dc=1+zeros(N,1);
x0.K=1+zeros(N,1);
x0.EbattV1=1.5+zeros(N,1);
options = optimoptions(prob.optimoptions,'Display','final');
% options = optimoptions(prob.optimoptions,'Algorithm','interior-point');
[values,fval,exitflag] = solve(prob,x0,'Options',options)

Solving problem using lsqnonlin.

Local minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the value of the optimality tolerance,
and constraints are satisfied to within the value of the constraint tolerance.

values = struct with fields:
    EbattV1: [10×1 double]
          K: [10×1 double]
        Pb1: [10×1 double]
      Pb1dc: [10×1 double]
     PgridV: [10×1 double]

fval = 58.5791

exitflag = 
    OptimalSolution

% Parse optmization results
if exitflag <= 0
    PgridV = zeros(N,1);
    Pb1 = zeros(N,1);
    Pb1dc = zeros(N,1);
    EbattV1 = zeros(N,1);
    K = zeros(N,1);
else
    PgridV = values.PgridV ;
    Pb1 =   values.Pb1;
    Pb1dc = values.Pb1dc
    EbattV1 = values.EbattV1;
    K = values.K    
end 

Pb1dc = 10×1

-0.9524 -0.9524 -0.3679 0.7251 -0.9524 1.0000 1.0000 1.0000 -0.9524 -0.0952

K = 10×1

0.9524 0.9524 0.3679 0.7251 0.9524 1.0000 1.0000 1.0000 0.9524 0.0952

[Jayachandra Malavatu]

i am using MATLAB 20222a, is this reason why it is used fimncon rather than lsqnonlin?

[Torsten]

Yes. Nonlinear constraints in lsqnonlin were introduced recently.