Optimization Formulation for "Complex" Constraints

Hello,

I am currently looking into solving an optimization problem using fmincon(). The problem is:

For this problem, fij's and tau_j's are the variables to be optimized. qobs_jk is a given, and q_jk is calculated as

where delta_t and I_ik are givens - the unknowns being optimized are tau and f.

My confusion stems from the following points/questions:

1. I am not sure how to treat the inputs tau and f into fmincon() because the objective function is not explicitly a function of tau or f.
2. How do I model the inequality constraint that contains a summation? Also, I need to satisfy this constraint for all i - can I do that in a single constraint?
3. Can the solver handle two (one of which is multi-dimensional) variable solutions? The optimal f_ij array should be of size (i x j) and the optimal tau_j array should be of size (j x 1).

I am fairly familiar with Matlab but brand new to using the optimization toolbox. Any help with formulating this problem using fmincon or guidance toward understanding how to formulate it correctly would be greatly appreciated.

Answer 1

[Matt J]

I am not sure how to treat the inputs tau and f into fmincon() because the objective function is not explicitly a function of tau or f.

Isn't it? In your formula for q_jk, all the tau's and f's appear on the right hand side. Anyway, the only thing that matters to fmincon is that you provide a function that will compute z for any given guess of the unknowns f and tau. It doesn't know or care what steps the function has to go through to reach z.

How do I model the inequality constraint that contains a summation?

That's what the Aineq,bineq arguments to fmincon are for.

Also, I need to satisfy this constraint for all i - can I do that in a single constraint?

No, you need n_i of them.

Can the solver handle two (one of which is multi-dimensional) variable solutions? The optimal f_ij array should be of size (i x j) and the optimal tau_j array should be of size (j x 1).

I would recommend, first of all, making a change of variables and replace exp(delta t/tau_j) with a single variable theta_j to simplify the computations and make them more linear. To answer your question though, f and tau (or f and theta) must be passed to your objective function bundled together in a single array unknowns. You can, of course, unpack them into separate variables within your objective function, if you wish.

To make it simpler to set this up, however, you can use some of the Optimization Toolbox tools for Problem-Based Optimization Setup together with the function prob2matrices() which must be downloaded. This will allow you to set up the optimization as below (where the blanks are for you to fill in):

f=optimvar('f',[ni,np],'LowerBound',0);

theta=optimvar('theta',np,'LowerBound',0,'UpperBound',1);

con.sumbound=sum(f,2)<=1;

sol0.f=____;  %Initial guesses
sol0.theta=___;

[p,idx]=prob2matrices({f,theta},'Constraints',con,'x0',sol0);

fun=@(unknowns) zfun(unknowns,idx, I, qobj);

  unknowns  = fmincon(fun,p.x0,p.Aineq,p.bineq,[],[],p.lb,p.ub);
  sol=x2sol(unknowns,idx);
  
f=sol.f;  %final solutions
theta=sol.theta;
  

function z=zfun(unknowns,idx, I, qobj)

  sol=x2sol(unknowns,idx);
  f=sol.f;
  theta=sol.theta;
  
  Q=f.'*I;
  
  q=nan(size(Q));
  q(:,1)=_____;
  for k=2:size(I,2)
      q(:,k)=q(:,k-1).*theta + (1-theta).*Q(:,k);
  end
      
  z=norm(qobj-q,'fro').^2;

end

[Corey Hoydic]

This is quite phenomenal, Matt. I will work on implementation and be sure to ask should I have any further questions or run into any problems.

[Matt J]

I hope it works out, but if it does, please do Accept-click the answer.

[Corey Hoydic]

Hello Matt,

I've got a couple questions so far:

1. I am trying to use your function prob2matrices [input as [p,idx]=prob2matrices({f,tau},'Constraints',con,'x0',sol0);], but I receive the error "'x0' is not a recognized parameter. For a list of valid name-value pair arguments, see the documentation for this function." Based on your documentation, I corrected the syntax to [p,idx]=prob2matrices({f,tau},'Constraints',con,'sol0',sol0);. However, I now receive the error

"Error using optim.internal.problemdef.ProblemImpl/mapSolution. MAPSOLUTION has been removed. Use VARINDEX instead." Using VARINDEX leads me to the error "Variable name must be a non-empty character vector or a string array." Do you have a suggestion to modify your function to make it functional with VARINDEX (I am using MATLAB R2021a)? Thanks!

Edit1: I see your question asked on the matter from Aug 2020, I will attempt to implement the solution but your insight is still welcomed.

Edit2: I think I fixed by replacing idx=mapSolution(prob,1:numel(allvars),[],[],[],[],[]); with idx=varindex(prob); (omitting the second argument). Can you confirm?

2. I am a bit apprehensive converting exp(-delta_t/tau) to theta, primarily because I see difficulty in establishing appropriate bounds (the bounds for tau are [0, inf]) coupled with the fact that delta_t may change between time steps. My formal NLP knowledge is limited, so I wonder, will omitting the variable substitution exp(-delta_t/tau)-->theta cause convergence problems, or is it merely for convenience's sake?

Really do appreciate your continued help. Thank you!

[Matt J]

I think I fixed by replacing idx=mapSolution(prob,1:numel(allvars),[],[],[],[],[]); with idx=varindex(prob); (omitting the second argument). Can you confirm?

Yes, that looks fine. [Sigh] I guess I'll have to update the FEX distribution.

I am a bit apprehensive converting exp(-delta_t/tau) to theta, primarily because I see difficulty in establishing appropriate bounds (the bounds for tau are [0, inf]) coupled with the fact that delta_t may change between time steps

As long as delta_t is positive, we know that exp(-delta_t/tau) will sweep through the interval (0,1) as tau sweep through [0, inf]. So, I don't think establishing bounds is the problem.

You can try both ways to see which works better. My main concern with writing the problem in terms of tau is that exp(-delta_t/tau) is highly flat as a function of tau near tau=0 (see below). Therefore, the solver may strain the limits of double float precision trying to compare different values of tau in that neighborhood.

fplot(@(tau) exp(-1./tau),[0,0.5])

[Corey Hoydic]

Thanks again, Matt. I have implemented the solution in terms of both tau and theta -- it would appear that both solutions give the same f's, but not the same tau's. Moreover, the solution for tau's seem to be sensitive to the initial guess for f and tau. This will require some further understanding of the "correct initial guess" on my part (though the f's are more important than the tau's for the work I am doing).

Thank you for all of your help, no more questions for now. Answer has been accepted.

[Matt J]

it would appear that both solutions give the same f's, but not the same tau's

But is the objective function approximately the same for both? If so, you simply have approximately non-unique solutions.

[Corey Hoydic]

The objective function is indeed approximately the same (for both theta and tau formulations AND for different initial guesses for f/tau/theta). The example problem I am using to debug this code is quite simple and idealized, so it will be interesting to see if there is non-uniqueness when applying this optimization to real data (I am working with oil well injection (I) and production (q) rates).