'fmincon', any way to improve the code?

Hello the community,

I was wondering if any one could have a look at my code and help me with the following questions.

1. Is it possible to keep "Theta" fixed at 0.02? With the current code, Theta keeps changing at each iteration.
2. I think my code is too long and too inefficient. I may have put irrelevant or unnecessary lines of code which make execution slow. Any suggestions on how to improve would also be appreciated!
3. Is it possible to pass Nc and Nw onto the function and run the code? at the moment, I have to give three initial values : Nc, Nw, and Lc. I want to give the first two only and Lc should come out of the calculations. just a long shot.

There are two functions, the first function should serve to allocate Nc and Nw values and the second function serves as an extra constraint for the first equation( function). I

clc
clear all
global sigma tauc tauw Theta c eta nu a a1 a2 A
format short
T=300;

a=0.6;       
a1=0.15;
a2=0.25;
tauw=0.5;        
tauc=1;        
Theta=0.02; % this should be constant

sigma=3;           
phic=3;             
phiw=6;            
gamma=1;

nu=10^3;            
eta=1;              
lambda=1;

%Initial values

Ac0=2;           
Aw0=10;

Sw0=4*10^3;
Sc0=35*10^3;
Scb=3.5*10^7;
Swb=1.5*10^7;
Ec=zeros(T,1); 
Ew=zeros(T,1);
X=zeros(T,1);
Ac=zeros(T+1,1); 
Aw=zeros(T+1,1);
Ac(1,1)=Ac0;
Aw(1,1)=Aw0;
Pc=zeros(T,1); 
Pw=zeros(T,1);
cc=zeros(T,1); 
cw=zeros(T,1);
Nc=zeros(T,1); 
Nw=zeros(T,1);
Sc=zeros(T+1,1); 
Sw=zeros(T+1,1);
Sc(1,1)=Sc0;    
Sw(1,1)=Sw0;
Lw=zeros(T,1);
Lc=zeros(T,1);
Xw=zeros(T,1);
Xc=zeros(T,1);
for t=1:T
A=[Ac(t),Aw(t)];

Xc(t)=Sc(t)/Ac(t)^lambda;
Xw(t)=Sw(t)/Aw(t)^lambda;
cc(t)=phic*(1+1/(Sc(t)+Xc(t)))^gamma;
cw(t)=phiw*(1+1/(Sc(t)+Xc(t)))^gamma;
c=[cc(t),cw(t)];

     x0=[500,500,0.5];

      UB=[ones(T,1); 100000*ones(T,1)]; % upper bound for the optimization
      LB=[zeros(T,1); zeros(T,1)]; % lower boud for the optimization
      options=optimset('TolFun',1e-11,'TolX',1e-9);
      out1=fmincon(@productivity_vL1,x0,[],[],[],[],LB,UB,@productivity_vL1c,options); 
      Nc(t,1)=out1(1);
      Nw(t,1)=out1(2);
      Lc(t,1)=out1(3);

      Ac(t+1)=(1+tauc*Theta)*Ac(t);
      Aw(t+1)=(1+tauw*Theta)*Aw(t);
  Lw(t)=1-Lc(t);
  Pc(t)=Lc(t)^(1/(1-sigma));
  Pw(t)=(1-Pc(t)^(1-sigma))^(1/(1-sigma));

Pc(t,1)=((Aw(t,1)/Ac(t,1))^(1-a1)*(cc(t,1)/cw(t,1))^a2)/(1+((Aw(t,1)/Ac(t,1))^(1-a1)*(cc(t,1)/cw(t,1))^a2)^(1-sigma))^(1/(1-sigma));
Pw(t,1)=(1-Pc(t,1)^(1-sigma))^(1/(1-sigma));
Lc(t,1)=Pc(t,1)^(1-sigma);
Lw(t,1)=1-Lc(t,1);

Ec(t)=(a2/cc(t)*Pc(t)^(1/(1-a1))*Lc(t)^((1-a)/(1-a1))*Ac(t+1))^((1-a1)/(1-a1*(1+a2)));
Ew(t)=(a2/cw(t)*Pw(t)^(1/(1-a1))*Lw(t)^((1-a)/(1-a1))*Aw(t+1))^((1-a1)/(1-a1*(1+a2)));

      %stock of energy
      Sc(t+1)=Sc(t)+Xc(t)-Ec(t);
      Sw(t+1)=Sw(t)+Xw(t)-Ew(t);

FUNCTION 2

function F=productivity_vL1(x)

    global A eta nu a a1 a2 c sigma tauc tauw Theta
    Nc=x(1);
    Nw=x(2);
    Lc=x(3);
    %
    Theta=nu*(Nc/((1+tauc)*A(1)))^eta;
    Theta=nu*(Nw/((1+tauc)*A(2)))^eta;
    Ac=(1+tauc*Theta)*A(1);
    Aw=(1+tauw*Theta)*A(2);

    Lw=1-Lc;
    Pc=Lc^(1/(1-sigma));
    Pw=(1-Pc^(1-sigma))^(1/(1-sigma));
    Ec=(a2/c(1)*Pc^(1/(1-a1))*Lc^((1-a)/(1-a1))*Ac)^((1-a1)/(1-a1*(1+a2)));
    Ew=(a2/c(2)*Pw^(1/(1-a1))*Lw^((1-a)/(1-a1))*Aw)^((1-a1)/(1-a1*(1+a2)));
    Lc=(Nc^(1-eta)/(nu*eta*a1*(1-a1)))^((1-a1)/(1-a))*(Pc*Ec^a2)^(1/(1-a));
    Lw=(Nw^(1-eta)/(nu*eta*a1*(1-a1)))^((1-a1)/(1-a))*(Pw*Ew^a2)^(1/(1-a));

    F=Lc+Lw-1;

FUNCTION2

function [f, ceq]=productivity_vL1c(x)

global A tauc tauw 

Nc=x(1);
Nw=x(2);
f=Nc/((1+tauc)*A(1))-Nw/((1+tauw)*A(2));
ceq=[];

Answer 1

[Matt J]

Is it possible to keep "Theta" fixed at 0.02? With the current code, Theta keeps changing at each iteration.

Of course, because you are deriving Theta from the unknowns, x(1) and x(2) which change with iterations. You have many constants in your code currently. Why not just set Theta to the desired constant and end of story?

I think my code is too long and too inefficient. I may have put irrelevant or unnecessary lines of code which make execution slow. Any suggestions on how to improve would also be appreciated!

One thing I see is that you are using a nonlinear constraint function productivity_vL1c to implement constraints that are in fact linear. Using the A,b inputs would be more efficient.

I want to give the first two only and Lc should come out of the calculations.

Again, why not just change your objective function code to calculate Lc the way you want? If you no longer want Lc to be an independent variable, you have control over the objective function, so you can derive Lc there any way you wish. Naturally, of course, you should then make x0 into a 2 variable vector instead of a 3 variable vector, etc...

[mrrox]

Thanks Matt,

Question 1. Here I want to keep Theta equal to 0.02 and let Nc and Nw vary. Therefore, if you noticed Theta is set 0.02 in the code, even then, Theta keeps changing. That's my problem for now.

Question 2. Yes, i realised it should not have been put in the 'nonlcon'.

Question 3. I have not yet found a way to do what you've just suggested. Still looking into it. Hope someone will help with rewriting the objective function if required. What I need to do is this, Nc and Nw should only be unknowns and the rest has to come out of computations. But the way the problem is formulated makes it difficult to avoid from making Lc an unknown. Lc and Lw should actually come out of the computation using the function.

Thank you for taking your time to help me. But would appreciate if you could help with the first question. Best Regards, R

[Matt J]

In what way does setting Theta=.02 constrain Nc and Nw? If no way, then why not just do this:

    function F=productivity_vL1(x)
        global A eta nu a a1 a2 c sigma tauc tauw Theta
        Nc=x(1);
        Nw=x(2);
        Lc=x(3);
        %
        Theta=.02;

        Ac=(1+tauc*Theta)*A(1);
        Aw=(1+tauw*Theta)*A(2);

        Lw=1-Lc;
        Pc=Lc^(1/(1-sigma));
        Pw=(1-Pc^(1-sigma))^(1/(1-sigma));
        Ec=(a2/c(1)*Pc^(1/(1-a1))*Lc^((1-a)/(1-a1))*Ac)^((1-a1)/(1-a1*(1+a2)));
        Ew=(a2/c(2)*Pw^(1/(1-a1))*Lw^((1-a)/(1-a1))*Aw)^((1-a1)/(1-a1*(1+a2)));
        Lc=(Nc^(1-eta)/(nu*eta*a1*(1-a1)))^((1-a1)/(1-a))*(Pc*Ec^a2)^(1/(1-a));
        Lw=(Nw^(1-eta)/(nu*eta*a1*(1-a1)))^((1-a1)/(1-a))*(Pw*Ew^a2)^(1/(1-a));
        F=Lc+Lw-1;

Otherwise, add an appropriate constraint on Nc and Nw.