Solving a DAE system: Error using decic function. index greater than 1 despite isLowIndexDAE = TRUE
I am trying to solve a system of DAEs that emerge in a constrained optimal control problem. I cannot find the consistent initial conditions. The error I get when looking for them using decic is that the index may be larger than 1, even though the isLowIndexDAE indicates that index is low. The error message is:
Error using decic>sls (line 170)
Index may be greater than one.
Error in decic (line 77)
[dy,dyp] = sls(res,dfdy,dfdyp,neq,free_y,free_yp);|
The variables psi1 and psi2 are the lagrange multipliers on two constraints: XB(t)>=0 and HW(t)=<1. Equations 6 and 7 are conseqiuently the complementary slackness conditions. I want to find the initial conditions in which both of these constraints bind, hence my guesses XB0=0 and HW0=1. Guesses for a few other variables follow. I tried relaxing the fixed guesses (and fix the relaxed ones) in the decic function in many different ways, but if I do that, then the only result is that I get a different error, for example:
Error using decic>sls (line 168)
Try freeing 7 fixed components.
Error in decic (line 77)
[dy,dyp] = sls(res,dfdy,dfdyp,neq,free_y,free_yp);
Error in ode_mathtoolbox (line 64)
[y0, yp0] = decic(F, 0, y0est, [0,1,1,1,0,1,1], yp0est, [0,0,1,1,0,0,0], opt)
My code is:
clear all
syms XB(t) HW(t) A(t) N(t) mu(t) psi1(t) psi2(t) …
rho nu n alpha alphahat chi lam phi delta s
eqn1 = N(t)/(A(t)*(1-s)*HW(t)*N(t)*alphahat – XB(t)) == (1-HW(t))*(1-nu)*(XB(t))^(-nu)+psi1(t);
eqn2 = (N(t)/(A(t)*(1-s)*HW(t)*N(t)*alphahat – XB(t)))*(A(t)*(1-s)*N(t)*alphahat) == XB(t)^(1-nu) – mu(t)*delta*lam*(A(t)^(phi))*((s*N(t))^lam)*(HW(t)^(lam-1)) + psi2(t);
eqn3 = diff(A(t), 1) == delta*(A(t)^phi)*(s*N(t)*HW(t))^lam;
eqn4 = diff(N(t), 1) == n*N(t);
eqn5 = diff(mu(t),1) == -(1-HW(t))*(1-nu)*(XB(t)^(-nu))*(1-s)*HW(t)*N(t)*alphahat + mu(t)*delta*phi*((A(t))^(phi-1))*(s*HW(t)*N)^lam + rho*mu(t);
eqn6 = XB(t)*psi1(t) == 0;
eqn7 = (1-HW(t))*psi2(t) == 0;
eqns = [eqn1 eqn2 eqn3 eqn4 eqn5 eqn6 eqn7];
vars = [XB(t), HW(t), A(t), N(t), mu(t), psi1(t), psi2(t)];
origVars = length(vars)
M = incidenceMatrix(eqns, vars)
isLowIndexDAE(eqns,vars)
[DAEs,DAEvars] = reduceDAEIndex(eqns,vars);
[DAEs,DAEvars] = reduceRedundancies(DAEs,DAEvars);
isLowIndexDAE(DAEs,DAEvars)
pDAEs = symvar(DAEs);
pDAEvars = symvar(DAEvars);
extraParams = setdiff(pDAEs, pDAEvars);
f = daeFunction(DAEs, DAEvars, alphahat, delta, lam, n, nu, phi, rho, s);
rho = 0.02;
nu = 0.9;
n = 0.01;
alpha = 0.33;
chi = 0.05;
lam = 0.65;
phi = 0.65;
delta = 0.057;
s = 0.2;
alphahat = alpha^(alpha/(1-alpha))-alpha^(1/(1-alpha));
F = @(t,Y,YP) f(t, Y, YP, alphahat, delta, lam, n, nu, phi, rho, s);
DAEvars
XB0 = 0;
HW0 = 1;
A0 = 1;
N0 = 1;
mu0 = 0;
y0est = [XB0; HW0; A0; N0; mu0; 1/(A0*(1-s)*N0*alphahat); 1+mu0*delta*lam*s*((A0)^phi)*(N0^lam)];
yp0est = [0,0,delta*(s^lam),n,0,0,0];
opt = odeset(‘RelTol’, 10.0^(-7), ‘AbsTol’ , 10.0^(-7));
[y0, yp0] = decic(F, 0, y0est, [0,0,0,1,0,0,0], yp0est, [0,0,0,0,0,0,0], opt)I am trying to solve a system of DAEs that emerge in a constrained optimal control problem. I cannot find the consistent initial conditions. The error I get when looking for them using decic is that the index may be larger than 1, even though the isLowIndexDAE indicates that index is low. The error message is:
Error using decic>sls (line 170)
Index may be greater than one.
Error in decic (line 77)
[dy,dyp] = sls(res,dfdy,dfdyp,neq,free_y,free_yp);|
The variables psi1 and psi2 are the lagrange multipliers on two constraints: XB(t)>=0 and HW(t)=<1. Equations 6 and 7 are conseqiuently the complementary slackness conditions. I want to find the initial conditions in which both of these constraints bind, hence my guesses XB0=0 and HW0=1. Guesses for a few other variables follow. I tried relaxing the fixed guesses (and fix the relaxed ones) in the decic function in many different ways, but if I do that, then the only result is that I get a different error, for example:
Error using decic>sls (line 168)
Try freeing 7 fixed components.
Error in decic (line 77)
[dy,dyp] = sls(res,dfdy,dfdyp,neq,free_y,free_yp);
Error in ode_mathtoolbox (line 64)
[y0, yp0] = decic(F, 0, y0est, [0,1,1,1,0,1,1], yp0est, [0,0,1,1,0,0,0], opt)
My code is:
clear all
syms XB(t) HW(t) A(t) N(t) mu(t) psi1(t) psi2(t) …
rho nu n alpha alphahat chi lam phi delta s
eqn1 = N(t)/(A(t)*(1-s)*HW(t)*N(t)*alphahat – XB(t)) == (1-HW(t))*(1-nu)*(XB(t))^(-nu)+psi1(t);
eqn2 = (N(t)/(A(t)*(1-s)*HW(t)*N(t)*alphahat – XB(t)))*(A(t)*(1-s)*N(t)*alphahat) == XB(t)^(1-nu) – mu(t)*delta*lam*(A(t)^(phi))*((s*N(t))^lam)*(HW(t)^(lam-1)) + psi2(t);
eqn3 = diff(A(t), 1) == delta*(A(t)^phi)*(s*N(t)*HW(t))^lam;
eqn4 = diff(N(t), 1) == n*N(t);
eqn5 = diff(mu(t),1) == -(1-HW(t))*(1-nu)*(XB(t)^(-nu))*(1-s)*HW(t)*N(t)*alphahat + mu(t)*delta*phi*((A(t))^(phi-1))*(s*HW(t)*N)^lam + rho*mu(t);
eqn6 = XB(t)*psi1(t) == 0;
eqn7 = (1-HW(t))*psi2(t) == 0;
eqns = [eqn1 eqn2 eqn3 eqn4 eqn5 eqn6 eqn7];
vars = [XB(t), HW(t), A(t), N(t), mu(t), psi1(t), psi2(t)];
origVars = length(vars)
M = incidenceMatrix(eqns, vars)
isLowIndexDAE(eqns,vars)
[DAEs,DAEvars] = reduceDAEIndex(eqns,vars);
[DAEs,DAEvars] = reduceRedundancies(DAEs,DAEvars);
isLowIndexDAE(DAEs,DAEvars)
pDAEs = symvar(DAEs);
pDAEvars = symvar(DAEvars);
extraParams = setdiff(pDAEs, pDAEvars);
f = daeFunction(DAEs, DAEvars, alphahat, delta, lam, n, nu, phi, rho, s);
rho = 0.02;
nu = 0.9;
n = 0.01;
alpha = 0.33;
chi = 0.05;
lam = 0.65;
phi = 0.65;
delta = 0.057;
s = 0.2;
alphahat = alpha^(alpha/(1-alpha))-alpha^(1/(1-alpha));
F = @(t,Y,YP) f(t, Y, YP, alphahat, delta, lam, n, nu, phi, rho, s);
DAEvars
XB0 = 0;
HW0 = 1;
A0 = 1;
N0 = 1;
mu0 = 0;
y0est = [XB0; HW0; A0; N0; mu0; 1/(A0*(1-s)*N0*alphahat); 1+mu0*delta*lam*s*((A0)^phi)*(N0^lam)];
yp0est = [0,0,delta*(s^lam),n,0,0,0];
opt = odeset(‘RelTol’, 10.0^(-7), ‘AbsTol’ , 10.0^(-7));
[y0, yp0] = decic(F, 0, y0est, [0,0,0,1,0,0,0], yp0est, [0,0,0,0,0,0,0], opt) I am trying to solve a system of DAEs that emerge in a constrained optimal control problem. I cannot find the consistent initial conditions. The error I get when looking for them using decic is that the index may be larger than 1, even though the isLowIndexDAE indicates that index is low. The error message is:
Error using decic>sls (line 170)
Index may be greater than one.
Error in decic (line 77)
[dy,dyp] = sls(res,dfdy,dfdyp,neq,free_y,free_yp);|
The variables psi1 and psi2 are the lagrange multipliers on two constraints: XB(t)>=0 and HW(t)=<1. Equations 6 and 7 are conseqiuently the complementary slackness conditions. I want to find the initial conditions in which both of these constraints bind, hence my guesses XB0=0 and HW0=1. Guesses for a few other variables follow. I tried relaxing the fixed guesses (and fix the relaxed ones) in the decic function in many different ways, but if I do that, then the only result is that I get a different error, for example:
Error using decic>sls (line 168)
Try freeing 7 fixed components.
Error in decic (line 77)
[dy,dyp] = sls(res,dfdy,dfdyp,neq,free_y,free_yp);
Error in ode_mathtoolbox (line 64)
[y0, yp0] = decic(F, 0, y0est, [0,1,1,1,0,1,1], yp0est, [0,0,1,1,0,0,0], opt)
My code is:
clear all
syms XB(t) HW(t) A(t) N(t) mu(t) psi1(t) psi2(t) …
rho nu n alpha alphahat chi lam phi delta s
eqn1 = N(t)/(A(t)*(1-s)*HW(t)*N(t)*alphahat – XB(t)) == (1-HW(t))*(1-nu)*(XB(t))^(-nu)+psi1(t);
eqn2 = (N(t)/(A(t)*(1-s)*HW(t)*N(t)*alphahat – XB(t)))*(A(t)*(1-s)*N(t)*alphahat) == XB(t)^(1-nu) – mu(t)*delta*lam*(A(t)^(phi))*((s*N(t))^lam)*(HW(t)^(lam-1)) + psi2(t);
eqn3 = diff(A(t), 1) == delta*(A(t)^phi)*(s*N(t)*HW(t))^lam;
eqn4 = diff(N(t), 1) == n*N(t);
eqn5 = diff(mu(t),1) == -(1-HW(t))*(1-nu)*(XB(t)^(-nu))*(1-s)*HW(t)*N(t)*alphahat + mu(t)*delta*phi*((A(t))^(phi-1))*(s*HW(t)*N)^lam + rho*mu(t);
eqn6 = XB(t)*psi1(t) == 0;
eqn7 = (1-HW(t))*psi2(t) == 0;
eqns = [eqn1 eqn2 eqn3 eqn4 eqn5 eqn6 eqn7];
vars = [XB(t), HW(t), A(t), N(t), mu(t), psi1(t), psi2(t)];
origVars = length(vars)
M = incidenceMatrix(eqns, vars)
isLowIndexDAE(eqns,vars)
[DAEs,DAEvars] = reduceDAEIndex(eqns,vars);
[DAEs,DAEvars] = reduceRedundancies(DAEs,DAEvars);
isLowIndexDAE(DAEs,DAEvars)
pDAEs = symvar(DAEs);
pDAEvars = symvar(DAEvars);
extraParams = setdiff(pDAEs, pDAEvars);
f = daeFunction(DAEs, DAEvars, alphahat, delta, lam, n, nu, phi, rho, s);
rho = 0.02;
nu = 0.9;
n = 0.01;
alpha = 0.33;
chi = 0.05;
lam = 0.65;
phi = 0.65;
delta = 0.057;
s = 0.2;
alphahat = alpha^(alpha/(1-alpha))-alpha^(1/(1-alpha));
F = @(t,Y,YP) f(t, Y, YP, alphahat, delta, lam, n, nu, phi, rho, s);
DAEvars
XB0 = 0;
HW0 = 1;
A0 = 1;
N0 = 1;
mu0 = 0;
y0est = [XB0; HW0; A0; N0; mu0; 1/(A0*(1-s)*N0*alphahat); 1+mu0*delta*lam*s*((A0)^phi)*(N0^lam)];
yp0est = [0,0,delta*(s^lam),n,0,0,0];
opt = odeset(‘RelTol’, 10.0^(-7), ‘AbsTol’ , 10.0^(-7));
[y0, yp0] = decic(F, 0, y0est, [0,0,0,1,0,0,0], yp0est, [0,0,0,0,0,0,0], opt) decic MATLAB Answers — New Questions
