Using minimization with variables in for loop

Hey guys, I just have a question about minimization and if you guys could please help me. I am pretty new to matlab and I know I have probably made some rookie mistakes.

I asked a question previously and I got no answer but I have then changed my code.

      function y = membrane(x,z)
      qpncal = 0;
      qpCO2cal = 0;
      t=1;
      N = 100;
      Am = 0.001010336;
      n = 10;
      Amn = Am/n;
      PCO2 = x;
      pf = 147256.4244;
      pp = 109138.049;
      xr(1) = 0.4718;
      yp(1) = 0.565;
      alpha = z;
      qp(1) = 0.000023078;
      qr(1) = 0.000046606;
      for j=1:N
          qCO2(j)= PCO2/t*Amn*((pf*xr(j))-(pp*yp(j)));
          qN2(j)= PCO2/t*Amn*(1/alpha)*((pf*(1-xr(j)))-(pp*(1-yp(j))));
          dqpn(j) = qCO2(j) + qN2(j);
          dypn(j) = qCO2(j)/(qCO2(j) + qN2(j));
          qr(j+1) = qr(j) + dqpn(j);
          xr(j+1) = (qr(j)*xr(j) + dqpn(j)*dypn(j))/qr(j + 1);
          qp(j+1) = qp(j) - dqpn(j);
          yp(j+1) = (qp(j)*yp(j) - dqpn(j)*dypn(j))/qp(j+1);
          qpncal = qpncal + dqpn(j);
          qpCO2cal = qpCO2cal + qCO2(j);
      end
        %0.58 is the value it has to converge to 
    y = qpCO2cal/qpncal-0.58
    end

So basically this function accepts two arguments, x being PCO2 and z being alpha. I want to minimize the function y and I am looking for the inputs of PCO2 and alpha that will do so. I tried running optimtool, putting fmincon, algorithm SQP, my objective function is @membrane, derivatives approximated by solver, start point [2 2], non-linear constraint function is @membrane, derivatives approximated by solver. I have read and seen people online doing the constraint function and optimization function separately and if somebody could please help me do so. Also when I run with the preceding settings, I get not enough input arguments. Could somebody please help?

Thank you.

Answer 1

[Walter Roberson]

Your nonlinear constraint function should not be the same as your objective function. Please read about nonlcon in the fmincon documentation. Note in particular that nonlcon is required to output two values; with you trying to use the single-output membrane function, you would get an error about not enough inputs.

[Ali]

Thanks a lot for taking your time and answering.

I now have two .m files (one for the function and one for the objective). Membrane is the objective and constraints is the constraints.

function constraints
      qpncal = 0;
      qpCO2cal = 0;
      t=1;
      N = 100;
      Am = 0.001010336;
      n = 10;
      Amn = Am/n;
      PCO2 = x;
      pf = 147256.4244;
      pp = 109138.049;
      xr(1) = 0.4718;
      yp(1) = 0.565;
      alpha = z;
      qp(1) = 0.000023078;
      qr(1) = 0.000046606;
      for j=1:N
          qCO2(j)= PCO2/t*Amn*((pf*xr(j))-(pp*yp(j)));
          qN2(j)= PCO2/t*Amn*(1/alpha)*((pf*(1-xr(j)))-(pp*(1-yp(j))));
          dqpn(j) = qCO2(j) + qN2(j);
          dypn(j) = qCO2(j)/(qCO2(j) + qN2(j));
          qr(j+1) = qr(j) + dqpn(j);
          xr(j+1) = (qr(j)*xr(j) + dqpn(j)*dypn(j))/qr(j + 1);
          qp(j+1) = qp(j) - dqpn(j);
          yp(j+1) = (qp(j)*yp(j) - dqpn(j)*dypn(j))/qp(j+1);
          qpncal = qpncal + dqpn(j);
          qpCO2cal = qpCO2cal + qCO2(j);
      end
end

function y = membrane(x,z)
    y = qpCO2cal/qpncal-0.58  
end

I open up optimtool, go to fmincon, interior point, objective function is @membrane, start point [1 1], and non-linear constraint function is @constraints. I now get

Optimization running.
Error running optimization.
In an assignment  A(I) = B, the number of elements in B and I must be the same.

I know you've mentioned that nonlcon is required to output two values and I think that is the error. Could you please clarify this a bit more? I want the two outputs to be the values of PCO2 and alpha that will minimize y.

Best regards, Ali

[Walter Roberson]

That isn't a constraint function! Not even close!

http://www.mathworks.com/help/toolbox/optim/ug/fmincon.html

The function that computes the nonlinear inequality constraints c(x)≤ 0 and the nonlinear equality constraints ceq(x) = 0. nonlcon accepts a vector x and returns the two vectors c and ceq. c is a vector that contains the nonlinear inequalities evaluated at x, and ceq is a vector that contains the nonlinear equalities evaluated at x. nonlcon should be specified as a function handle to a file or to an anonymous function, such as mycon:
x = fmincon(@myfun,x0,A,b,Aeq,beq,lb,ub,@mycon)
where mycon is a MATLAB function such as
function [c,ceq] = mycon(x)
c = ...     % Compute nonlinear inequalities at x.
ceq = ...   % Compute nonlinear equalities at x.

I do not, however, see any constraints described in your problem description. What value limits are there? For example, does alpha need to be non-negative? It is not obvious to me that you have any constraints at all, so it is not obvious to me that you should be using fmincon()

You also have the problem that you are trying to minimize over two separate variables by passing the two variables in separately. The function to be evaluated, the one whose handle is the first argument to fmincon, must accept only a single variable (that can be a vector.) If you want to minimize over the two scalars x and z, then make the first argument of the function a vector and break apart the vector at the top of the function:

function y = membrane(x)
  PC02 = x(1);
  z = x(2);
  qpncal = 0;
  qpCO2cal = 0;
  t=1;
  %[etc]
  y = qpCO2cal/qpncal-0.58;
end

[Ali]

Okay, I had that initially (x(1) and x(2)).

Alpha has to be >0 and PCO2 has to be >0, if we can include that in the constraint function.

I have then changed my code, I want to use fsolve to solve the following two equations simultaneously.

function y = membrane(x)
qpncal = 0;
qpCO2cal = 0;
t=1;
N = 100;
Am = 0.001010336;
n = 10;
Amn = Am/n;
PCO2 = x(1);
pf = 147256.4244;
pp = 109138.049;
xr(1) = 0.4718;
yp(1) = 0.565;
alpha = x(2);
qp(1) = 0.000023078;
qr(1) = 0.000046606;
for j=1:N
    qCO2(j)= PCO2/t*Amn*((pf*xr(j))-(pp*yp(j)));
    qN2(j)= PCO2/t*Amn*(1/alpha)*((pf*(1-xr(j)))-(pp*(1-yp(j))));
    dqpn(j) = qCO2(j) + qN2(j);
    dypn(j) = qCO2(j)/(qCO2(j) + qN2(j));
    qr(j+1) = qr(j) + dqpn(j);
    xr(j+1) = (qr(j)*xr(j) + dqpn(j)*dypn(j))/qr(j + 1);
    qp(j+1) = qp(j) - dqpn(j);
    yp(j+1) = (qp(j)*yp(j) - dqpn(j)*dypn(j))/qp(j+1);
    qpncal = qpncal + dqpn(j);
    qpCO2cal = qpCO2cal + qCO2(j);
end
  %0.58 is the value it has to converge to 
    y = [ qpncal - 2;
          qpCO2cal/qpncal-0.58;];
end

I make a vector x0 = [1 1] in the command window. Then I call x = fsolve(@membrane, x0) and I get the following error.

In an assignment  A(I) = B, the number of elements in B and I must be the same.

Error in membrane (line 53)
   alfa(j) = alfa(j-1) + del;

Error in fsolve (line 241)
            fuser = feval(funfcn{3},x,varargin{:});

Caused by:
    Failure in initial user-supplied objective function evaluation. FSOLVE cannot
    continue.

I specifically don't understand the first statement, In an assignment.... etc...

Thanks a lot for the help.

Answer 2

[Ali]

I get some promising results now. I ran the function as x = fsolve(@membrane, [2e-6 1]) and I get

Equation solved, fsolve stalled.

fsolve stopped because the relative size of the current step is less than the
default value of the step size tolerance squared and the vector of function values
is near zero as measured by the default value of the function tolerance.

    <stopping criteria details>
x =

    0.0000    1.0000

How can you increase decimal accuracy?

[Ali]

I read online and you can do >>format long and this works. Thanks a lot Walter for your input and I was wondering if there is any way you can increase the accuracy of the answer. It's okay if it takes a longer time to iterate but I want accuracy.