Multivariate Nonlinear Problem With vector & matrix unknown inputs

Hello,

Problem description:

I have oil flow rate observed (qo_obs) and I have a model to calculate oil flow rate (qo_cal).

The model to find qo_cal is somewhat complex, and it invloved three unknowns: f, tau, and V.

f is a 9x16 matrix, tau is a 9x1 vector, V is a 9x1 vector.

each of the unkowns above has its own constraints, such as:

Sum of every column of f = 1, (How do I even write this constraint?)

any tau >0

0< V<C, where C is some constant I know.

Now I want to find the values of these unkowns that will yield minimum of (qo_obs-qo_cal)^2

What solver to use? and how?

Thanks,

Answer 1

[Matt J]

It's very easy to set up all the variables and constraints with the problem-based approach.

f=optimvar('f',[9,16]);
tau=optimvar('tau',[9,1],'LowerBound',zeros(9,1));
V=optimvar('V',[9,1],'LowerBound',0,'UpperBound',C);

prob=optimproblem('ObjectiveSense','minimize');
prob.Constraints.fconstraint=sum(f)==1;

prob.Objective=_____________

Once you've defined your objective, just do

sol=solve(prob)

[Ahmed Ghamdi]

Thanks Matt,

Can you explain please:

prob.Constraints.fconstraint=sum(f)==1;

The above line does it specify that each column in f must equal 1? Because this is the constraint I want to impose.

prob.Objective=_____________

What should I make my objective function? the square difference?

How do I check the values for the optimized unknowns? whenI click on them in workspace it tells me that they became OptimizationVariables, but did not show me any values.

finally, the optimization gave the following message 'Solving problem using linprog'. This is a linear regression I believe, is there a way to make it nonlinear?

Thanks,

[Matt J]

Well, first read the link I gave you on the problem-based approach. That will answer a lot of these questions and make follow-up conversation easier.

[Ahmed Ghamdi]

I have read most of the problem-based approach link, I only have one issue remaining, is that how do I do any operation with the OptimizationVariables once they optimization completed. when I tried to do for example f(1,1)+f(1,2) is shows me

ans =

Linear OptimizationExpression

f(1, 1) + f(1, 2)

but what are the values? as in double values?

at least I've gotten to this point in my optimization problem, which is way ahead of when I began early today, so thank you, I will try to figure it out.

[Matt J]

I suggest you show the code you have used to generate f. The output of the optimization should be doubles.

[Ahmed Ghamdi]

Hey Matt, thanks for following up,

I have found out how to get the values of f, I need to type

sol.f

and if I want a particular value of f, say f(1,1) I type

sol.f(1,1)

But my problem is that now I try to find values of the other optimized values, tau and Vp, when I type sol.tau or sol.Vp I get

>> sol.tau
Reference to non-existent field 'tau'.

>> sol.Vp
Reference to non-existent field 'Vp'.

Below is my code

%% Assume values for f & tau
f = zeros(9,16);
tau = zeros(9,1);
for i=1:9
    tau(i,1) = 10;
    for j=1:16
    f(i,j) = 0.1111;
    end
end
%% Assume values for Vp
Vp = zeros(9,1);
for j=1:9
        Vp(j) = 10^5;
end

%% Start calculating oil rate
for j=1:9
    for k=2:120
        krw(k,j) = krw_endpoint*((Sw(k,j)-Swr)/(1-Swr-Sor))^n1;
        kro(k,j) = kro_endpoint*((1-Sw(k,j)-Sor)/(1-Swr-Sor))^n2;
        fw(k,j) = (krw(k,j)/uw)/((krw(k,j)/uw)+(kro(k,j)/uo));
        fo(k,j) = 1-fw(k,j);
        Mt(k,j) = krw(k,j)/uw+kro(k,j)/uo;
        for i=1:16
            press_sup(i) = f(j,i)*All_Inj_qw(k,i); %here i use f
        end
        summation = sum(press_sup);
        qt(k,j) = qt(k-1,j)*exp((-dt(k))/(tau(j)/Mt(k,j)))+(1-exp((-dt(k))/(tau(j)/Mt(k,j))))*summation; %here I use tau
        qo(k,j) = qt(k,j)*fo(k,j);
        So_avg(k,j) = So_avg(k-1,j)-dt(k)/Vp(j)*(((So_avg(k-1,j)*(cf+co))/ct)*(summation-qt(k,j))-qo(k,j)); %here I use Vp
        Q(k) = (Q(k-1)+summation)/Vp(j);
        So_outlet(k,j) = So_avg(k-1,j)+Q(k)*(1-fw(k,j));
        Sw(k+1,j) = 1-So_outlet(k,j);
    end
end

%% Calculate Error

for j=1:9
    for k=1:120
        Error(k,j) = (qo(k,j)-All_Prod_qo(k,j))^2;
    end
end
Total_Error = sum(Error);

%% Optimization Problem for finding optimized f,tau, and Vp

prob=optimproblem('ObjectiveSense','minimize');
f=optimvar('f',[9,16],'Type','continuous','LowerBound',0,'UpperBound',1);
tau=optimvar('tau',[9,1],'Type','continuous','LowerBound',zeros(9,1));
Vp=optimvar('Vp',[9,1],'Type','continuous','LowerBound',0,'UpperBound',1.6895e+07);
prob.Constraints.fconstraint=sum(f)==1;

prob.Objective = Total_Error;

sol = solve(prob);

Also, when I check my problem by typing

show(prob)

this what comes out

>> show(prob)

  OptimizationProblem : 

	Solve for:
       f

	minimize :
       1703960721.818
       

why isn't MATLAB solving for tau and Vp, why only f?

[Matt J]

It doesn't matter - the solution you have obtained is meaningless, because you have set prob.Objective to a constant number, 1703960721.818. You have used the f optimvar to define the constraint fconstraint, but you have not used any of the optimvar variables that you have created to define prob.Objective.

[Ahmed Ghamdi]

I have set prob.Objective to Total_Error which is the sum of square errors between qo_calc and qo__obs. And qo_calc is an expression linked to the optimvars variables I have.

[Matt J]

No, the lines in the code that create your optimization variables are these,

f=optimvar('f',[9,16],'Type','continuous','LowerBound',0,'UpperBound',1);
tau=optimvar('tau',[9,1],'Type','continuous','LowerBound',zeros(9,1));
Vp=optimvar('Vp',[9,1],'Type','continuous','LowerBound',0,'UpperBound',1.6895e+07);

You have not used these variables to generate Total_Error. You have used dummy numeric values defined above these lines.

[Matt J]

Ahmed's comment moved here:

Excellent, I see what you mean,

But now if I define and use the optimvar variables for calculating Total_Error, one of them (tau) is in an exponential function which gives me error. Any way around this?

By the way, I really appreciate your help,

[Matt J]

It appears that the form of your total error function is outside the scope of the problem based solver. However, you can convert the prob object to a form that fmincon can handle. This will look like below. You will have to implement your objective as an ordinary function,

myprob=prob2struct(prob);

myprob.f=@objective;

x=fmincon(myprob)

function TotalError= objective(x)

  f=x(1:9*16);
  tau=x(9*16+1:9*16+9)
  V=x(9*16+10:9*16+18)

  ...


end

[Matt J]

Ahmed commented

Thanks, i will try this when I get home.

By the way, how did you get all this MATLAB knowledge? What references do you suggest that will help me learn optimization in MATLAB?

[Matt J]

These are a good starting point

https://www.mathworks.com/help/optim/nonlinear-programming.html

https://www.mathworks.com/help/optim/ug/optimization-toolbox-tutorial.html

[Ahmed Ghamdi]

Thanks Matt,

my code is working fine, but fmincon takes 7 hours to run. Is this normal? Also, are there more advanced algorathims better than fmincon?

[Matt J]

There is no normal. It depends in part on how much time it takes to execute your Objective function code. I see in your code that you have done no vectorization at all - everything is done with for-loops. It is plausible to me that that would make execution time very long.