search

do-mpc / do-mpc

Model predictive control python toolbox
https://www.do-mpc.com/
GNU Lesser General Public License v3.0
1k stars 181 forks source link

IndexError: index 0 is out of bounds for axis 0 with size 0 #30

Closed biswa007753 closed 4 years ago

biswa007753 commented 4 years ago

Trying to implement the code in a hybrid electric vehicle but getting the following error in Spyder python 3.7

File "C:\ProgramData\Anaconda3\lib\site-packages\spyder_kernels\customize\spydercustomize.py", line 668, in runfile execfile(filename, namespace)

File "C:\ProgramData\Anaconda3\lib\site-packages\spyder_kernels\customize\spydercustomize.py", line 108, in execfile exec(compile(f.read(), filename, 'exec'), namespace)

File "C:/Users/biswajitt/Downloads/do-mpc-master/do-mpc-master/examples/batch_reactor/do-mpc.py", line 65, in configuration_1.setup_solver()

File "../../code\core_do_mpc.py", line 190, in setup_solver nlp_dict_out = setup_nlp.setup_nlp(self.model, self.optimizer)

File "../../code\setup_nlp.py", line 516, in setup_nlp P_ksb = P[b + branch_offset[k][s]]

IndexError: index 0 is out of bounds for axis 0 with size 0

alexandru-tc commented 4 years ago

Hello and thanks for showing interest in do-mpc!

It sounds like you are not defining the "uncertain parameters" in the model of your MPC. By default, we consider at least 1 parameter to be defined in the model template, and the values to be defined in the optimizer template. Please do this, even if you do not use the parameter in the model.Otherwise you will get errors during the setup of the NLP. We have implemented a very general NLP definition, that's the reason behind this :)

I suggest to look at the CSTR example and to use something similar to that example, although there we define 2 uncertain parameters. Hope this helps. Let me know if you have further questions.

Regards,

On Mon, Dec 2, 2019 at 9:48 AM biswa007753 notifications@github.com wrote:

Trying to implement the code in a hybrid electric vehicle but getting the following error in Spyder python 3.7

File "C:\ProgramData\Anaconda3\lib\site-packages\spyder_kernels\customize\spydercustomize.py", line 668, in runfile execfile(filename, namespace)

File "C:\ProgramData\Anaconda3\lib\site-packages\spyder_kernels\customize\spydercustomize.py", line 108, in execfile exec(compile(f.read(), filename, 'exec'), namespace)

File "C:/Users/biswajitt/Downloads/do-mpc-master/do-mpc-master/examples/batch_reactor/do-mpc.py", line 65, in configuration_1.setup_solver()

File "../../code\core_do_mpc.py", line 190, in setup_solver nlp_dict_out = setup_nlp.setup_nlp(self.model, self.optimizer)

File "../../code\setup_nlp.py", line 516, in setup_nlp P_ksb = P[b + branch_offset[k][s]]

IndexError: index 0 is out of bounds for axis 0 with size 0

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/do-mpc/do-mpc/issues/30?email_source=notifications&email_token=ACNIPGGC234Y56DHPTAP7ITQWTDV3A5CNFSM4JTSJFH2YY3PNVWWK3TUL52HS4DFUVEXG43VMWVGG33NNVSW45C7NFSM4H5GCYUQ, or unsubscribe https://github.com/notifications/unsubscribe-auth/ACNIPGDFFS4NC62LERTTDLDQWTDV3ANCNFSM4JTSJFHQ .

biswa007753 commented 4 years ago

Thamks for the help

biswa007753 commented 4 years ago

Hi I am unable to find a cost function to the following code as they are not converging if you could help

#

This file is part of do-mpc

#

do-mpc: An environment for the easy, modular and efficient implementation of

robust nonlinear model predictive control

#

Copyright (c) 2014-2018 Sergio Lucia, Alexandru Tatulea-Codrean

TU Dortmund. All rights reserved

#

do-mpc is free software: you can redistribute it and/or modify

it under the terms of the GNU Lesser General Public License as

published by the Free Software Foundation, either version 3

of the License, or (at your option) any later version.

#

do-mpc is distributed in the hope that it will be useful,

but WITHOUT ANY WARRANTY; without even the implied warranty of

MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the

GNU Lesser General Public License for more details.

#

You should have received a copy of the GNU General Public License

along with do-mpc. If not, see http://www.gnu.org/licenses/.

# from casadi import * import numpy as NP import core_do_mpc

def model():

"""
--------------------------------------------------------------------------
template_model: define the non-uncertain parameters
--------------------------------------------------------------------------
"""
M_value = M = 0.5;
m_value = m = 0.2;
b_value = b = 0.1;
I_value = I = 0.006;
g_value = g = 9.8;
l_value=  l = 0.3;
Vspeed= 20
"""
--------------------------------------------------------------------------
template_model: define uncertain parameters, states and controls as symbols
--------------------------------------------------------------------------
"""
# Define the uncertain/varying parameters as CasADi symbols
alpha      = SX.sym("alpha")         # can be used for implementing uncertainty

tv_param_1 = SX.sym("tv_param_1")    # can be used to implement setpoint change
x =SX.sym("x")
v =SX.sym("v")

# Define the differential states as CasADi symbols
dd = SX.sym("dd",4)
dd[0]  = -(I_value+m_value*(l_value**2))* b_value/I_value*(M_value+m_value)+M_value*m_value*l_value**2
dd[1]  = m_value**2*g_value*l_value**2/I_value*( M_value+m_value)+M_value*m_value*l_value**2
dd[2] = -m_value*l_value*b_value/ I_value*(M_value+m_value)+M_value*m_value*l_value
dd[3] = m_value**2*g_value*l_value*(M_value+m_value)/I_value*(M_value+m_value)+M_value*m_value*l_value**2          # angular velocity of the mass
# Define the algebraic states as CasADi symbols
#y      = SX.sym("y")          # the y coordinate is considered algebraic

# Define the control inputs as CasADi symbols
F   = SX.sym("F")             # control force applied to the lever

"""
--------------------------------------------------------------------------
template_model: define algebraic and differential equations
--------------------------------------------------------------------------
"""
# Define the differential equations
a = dd[0]    = SX.sym("dd[0]")               # x position of the mass
b =dd[1]       = SX.sym("dd[1]")                  # y position of the mass
c = dd[2]     = SX.sym("dd[2]")              # linear velocity of the mass
theta = SX.sym("theta")          # angle of the metal rod
e=dd[3] = SX.sym("dd[3]") 
# Define the algebraic equations
            # the coordinates must fulfil the constraint of the lever arm length

# Concatenate differential states, algebraic states, control inputs and right-hand-sides

_x = vertcat(a,b,c,e)

_z = []                                # toggle if there are no AE in your model
#_z = vertcat(y)

_u = vertcat(F)

_p = vertcat(alpha)

_xdot = vertcat(dd)

_zdot = []                          # toggle if there are no AE in your model
#_zdot = vertcat(dtrajectory)

_tv_p = vertcat(tv_param_1)
"""
--------------------------------------------------------------------------
template_model: initial condition and constraints
--------------------------------------------------------------------------
"""
# Initial conditions for Differential States
x_init  = 2.5
v_init  = 0.0
t_init  = 1.5
o_init  = 0.0

#Initial conditions for Algebraic States
y_init  = sqrt(3.0)/2.0

x0 = NP.array([x_init, v_init, t_init, o_init])
z0 = NP.array([])
# Bounds on the states. Use "inf" for unconstrained states
x_lb    =  -10.0;       x_ub  = 10.0
v_lb    = -10.0;       v_ub  = 10.0
t_lb    =   -100*pi;       t_ub  = 100*pi
o_lb    = -10.0;       o_ub  = 10.0

x_lb = NP.array([x_lb, v_lb, t_lb, o_lb])
x_ub = NP.array([x_ub, v_ub, t_ub, o_ub])
z_lb = NP.array([])
z_ub = NP.array([])
# Bounds on the control inputs. Use "inf" for unconstrained inputs
F_lb    = -250.0;       F_ub = 250.00 ;     F_init = 0.0    ;

u_lb=NP.array([F_lb])
u_ub=NP.array([F_ub])
u0 = NP.array([F_init])

# Scaling factors for the states and control inputs. Important if the system is ill-conditioned
x_scaling = NP.array([1.0, 1.0, 1.0, 1.0])
z_scaling = NP.array([])
u_scaling = NP.array([1.0])

# Other possibly nonlinear constraints in the form cons(x,u,p) <= cons_ub
# Define the expresion of the constraint (leave it empty if not necessary)
cons = vertcat([])
# Define the lower and upper bounds of the constraint (leave it empty if not necessary)
cons_ub = NP.array([])

# Activate if the nonlinear constraints should be implemented as soft constraints
soft_constraint = 0
# Penalty term to add in the cost function for the constraints (it should be the same size as cons)
penalty_term_cons = NP.array([])
# Maximum violation for the constraints
maximum_violation = NP.array([0])

# Define the terminal constraint (leave it empty if not necessary)
cons_terminal = vertcat([])
# Define the lower and upper bounds of the constraint (leave it empty if not necessary)
cons_terminal_lb = NP.array([])
cons_terminal_ub = NP.array([])

"""
--------------------------------------------------------------------------
template_model: cost function
--------------------------------------------------------------------------
"""
# Define the cost function
# Lagrange term
##-----Distance to vertical position---------------------
lterm =  (v-Vspeed)**2 + (theta+pi/2)**2
# Mayer term
mterm =  0
# Penalty term for the control movements
rterm = 0.0001*NP.array([1.0])

"""
--------------------------------------------------------------------------
template_model: pass information (not necessary to edit)
--------------------------------------------------------------------------
 """
model_dict = {'x':_x,'u': _u, 'rhs':_xdot,'p': _p, 'z':_z,'x0': x0,'x_lb': x_lb,'x_ub': x_ub, 'u0':u0, 'u_lb':u_lb, 'u_ub':u_ub, 'x_scaling':x_scaling, 'u_scaling':u_scaling, 'cons':cons,
"cons_ub": cons_ub, 'cons_terminal':cons_terminal, 'cons_terminal_lb': cons_terminal_lb,'tv_p':_tv_p, 'cons_terminal_ub':cons_terminal_ub, 'soft_constraint': soft_constraint, 'penalty_term_cons': penalty_term_cons, 'maximum_violation': maximum_violation, 'mterm': mterm,'lterm':lterm, 'rterm':rterm}

model = core_do_mpc.model(model_dict)

return model
alexandru-tc commented 4 years ago

It is quite hard to see what you are trying to achieve with the implementation. For an inverted cart with pendulum, the typical target state is the vertical position where the controller should stabilize the mass. If there is no external disturbance, the system should stop after achieving that state. If you're trying to target a different position (x,y) and/or angle theta and/or linear velocities, you should check if that point is a feasible one. Maybe if you provide more information, I could help you further. Other than that, the implementation should be okay. Just have a look at the original pendulum implementation and start from there.

Cheers! Alex

On Tue, Dec 3, 2019 at 10:20 AM biswa007753 notifications@github.com wrote:

Hi, I am unable to find a required cost function as they are unable to converge provided in the below cart pole inverted pendulum system

from casadi import * import numpy as NP import core_do_mpc

def model():

"""

template_model: define the non-uncertain parameters

""" m = 1.0 # mass of the ball [kg] l = 1.0 # length of lever arm [m] h = 0.582 # height of rod connection point [m] g = 9.81 # grav acceleration [m/s^2] Kx = 1 Txx = 1 Ky = 0.29 Txy = 0.37 theta=pi/10 theta1=pi/10 theta2=pi/10 """

template_model: define uncertain parameters, states and controls as symbols

"""

uncertain/varying paDefine therameters as CasADi symbols

     # can be used for implementing uncertainty

gamma = SX.sym("gamma") tv_param_1 = SX.sym("tv_param_1") beta = SX.sym("beta")

tv_param_1 = SX.sym("tv_param_1") # can be used to implement setpoint change

Define the differential states as CasADi symbols

alpha =SX.sym("alpha") speedx =SX.sym("speedx") positionx =SX.sym("positionx") speedy =SX.sym("speedy") positiony =SX.sym("positiony") theta1=SX.sym("theta1") theta2=SX.sym("theta2")

         # angular velocity of the mass

Define the algebraic states as CasADi symbols

     #y= SX.sym("y")          # the y coordinate is considered algebraic

Define the control inputs as CasADi symbols

F = SX.sym("F") # control force applied to the lever

"""

template_model: define algebraic and differential equations

"""

Define the differential equations

dd = SX.sym("dd",6) dd[0] = 0.714gsin(alpha) dd[1] = Kx/Txx*(theta1-alpha) dd[2] = speedx

dd[3] = 0.714gsin(beta) dd[4] = Kx/Txx*(theta2-beta) dd[5] = speedy

Concatenate differential states, algebraic states, control inputs and right-hand-sides

_x = vertcat(speedx,alpha,positionx,speedy,beta,positiony) _z = [] # toggle if there are no AE in your model

_z = vertcat(y)

_u = vertcat(theta1,theta2)

_p = vertcat(gamma)

_xdot = vertcat(dd)

_zdot = [] # toggle if there are no AE in your model

_zdot = vertcat(dtrajectory)

_tv_p = vertcat(tv_param_1) """

template_model: initial condition and constraints

"""

Initial conditions for Differential States

speedx_init = 0.1 alpha_init = pi/10 positionx_init = 2.0 speedy_init = 0.2 beta_init = pi/10 positiony_init = 2.0

x0 = NP.array([speedx_init,alpha_init,positionx_init,speedy_init,beta_init,positiony_init]);

Bounds on the states. Use "inf" for unconstrained states

speedx_lb = -1.0; speedx_ub = 1.0; alpha_lb = -pi/18.0; alpha_ub = pi/18.0; positionx_lb = -10.0; positionx_ub = 10.0;

speedy_lb = -1.0; speedy_ub = 1.0; beta_lb = -pi/18; beta_ub = pi/18; positiony_lb = -10.0; positiony_ub = 10.0;

s_lb = NP.array([speedx_lb, alpha_lb,positionx_lb,speedy_lb,beta_lb,positiony_lb]) s_ub = NP.array([speedx_ub, alpha_ub,positionx_ub,speedy_ub,beta_ub,positiony_ub])

Bounds on the control inputs. Use "inf" for unconstrained inputs

theta1_lb = -pi/20; theta1_ub = pi/20 ; theta1_init = pi/30 ; theta2_lb = -pi/20; theta2_ub = pi/20 ; theta2_init = pi/30 ;

u_lb=NP.array([theta1_lb,theta2_lb]) u_ub=NP.array([theta1_ub,theta2_ub]) u0 = NP.array([theta1_init,theta2_init])

Scaling factors for the states and control inputs. Important if the system is ill-conditioned

x_scaling = NP.array([1.0, 1.0,1.0,1.0,1.0,1.0]) z_scaling = NP.array([]) u_scaling = NP.array([1.0,1.0])

Other possibly nonlinear constraints in the form cons(x,u,p) <= cons_ub

Define the expresion of the constraint (leave it empty if not necessary)

cons = vertcat([])

Define the lower and upper bounds of the constraint (leave it empty if not necessary)

cons_ub = NP.array([])

Activate if the nonlinear constraints should be implemented as soft constraints

soft_constraint = 0

Penalty term to add in the cost function for the constraints (it should be the same size as cons)

penalty_term_cons = NP.array([])

Maximum violation for the constraints

maximum_violation = NP.array([0])

Define the terminal constraint (leave it empty if not necessary)

cons_terminal = vertcat([])

Define the lower and upper bounds of the constraint (leave it empty if not necessary)

cons_terminal_lb = NP.array([]) cons_terminal_ub = NP.array([])

"""

template_model: cost function

"""

Define the cost function

Lagrange term

-----Distance to vertical position---------------------

lterm = (speedx)2+(speedy)2 +(positionx-3)2+(positiony-4)2

Mayer term

mterm = 0

Penalty term for the control movements

rterm = 0.0001*NP.array([1.0])

"""

template_model: pass information (not necessary to edit)

""" model_dict = {'x':_x,'u': _u, 'rhs':_xdot,'p': _p, 'z':_z,'x0': x0,'x_lb': s_lb,'x_ub': s_ub, 'u0':u0, 'u_lb':u_lb, 'u_ub':u_ub, 'x_scaling':x_scaling, 'u_scaling':u_scaling, 'cons':cons, "cons_ub": cons_ub, 'cons_terminal':cons_terminal, 'cons_terminal_lb': cons_terminal_lb,'tv_p':_tv_p, 'cons_terminal_ub':cons_terminal_ub, 'soft_constraint': soft_constraint, 'penalty_term_cons': penalty_term_cons, 'maximum_violation': maximum_violation, 'mterm': mterm,'lterm':lterm, 'rterm':rterm}

model = core_do_mpc.model(model_dict)

return model

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/do-mpc/do-mpc/issues/30?email_source=notifications&email_token=ACNIPGANXWZPO4HBFFEJINTQWYQEHA5CNFSM4JTSJFH2YY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEFYVEAI#issuecomment-561074689, or unsubscribe https://github.com/notifications/unsubscribe-auth/ACNIPGDXEVQG4A6HZ34JFDTQWYQEHANCNFSM4JTSJFHQ .

biswa007753 commented 4 years ago

Hi Alex, the equation implemented is given in the following link http://ctms.engin.umich.edu/CTMS/index.php?example=InvertedPendulum&section=SystemModeling if you can implement it will be nice. I have implemented it in the above code My major concern is the cost function

biswa007753 commented 4 years ago

Hi Alex, another error which i am unable to resolve

File "C:\ProgramData\Anaconda3\lib\site-packages\casadi\casadi.py", line 14336, in integrator return _casadi.integrator(*args)

RuntimeError: Error in Function::Function for 'dae' [SXFunction] at .../casadi/core/function.cpp:219: .../casadi/core/function_internal.cpp:124: Error calling SXFunction::init for 'dae': .../casadi/core/x_function.hpp:238: Xfunction input arguments must be purely symbolic. Argument 1(x) is not symbolic.

The code implemented for the above error is below and the code is for a DC Motor implementation

from casadi import * import numpy as NP import core_do_mpc

def model():

"""
--------------------------------------------------------------------------
template_model: define the non-uncertain parameters
--------------------------------------------------------------------------
"""
J = 0.01;
b = 0.1;
K = 0.01;
R = 1;
L = 0.5;
"""
--------------------------------------------------------------------------
template_model: define uncertain parameters, states and controls as symbols
--------------------------------------------------------------------------
"""
#  uncertain/varying paDefine therameters as CasADi symbols
         # can be used for implementing uncertainty
speed  =s     = SX.sym("s")
tv_param_1 = SX.sym("tv_param_1")    # can be used to implement setpoint change

# Define the differential states as CasADi symbols

             # angular velocity of the mass
# Define the algebraic states as CasADi symbols
         #y= SX.sym("y")          # the y coordinate is considered algebraic

# Define the control inputs as CasADi symbols
v   = SX.sym("v")             # control force applied to the lever

"""
--------------------------------------------------------------------------
template_model: define algebraic and differential equations
--------------------------------------------------------------------------
"""
# Define the differential equations
dd = SX.sym("dd",2)
theta=dd[0] = -(b/J)+K/J
i= dd[1] = -(K/L)-(R/L)

# Concatenate differential states, algebraic states, control inputs and right-hand-sides

_x = vertcat(theta,i)
_z = []                                # toggle if there are no AE in your model
#_z = vertcat(y)

_u = vertcat(v)

_p = vertcat()

_xdot = vertcat(dd)

_zdot = []                          # toggle if there are no AE in your model
#_zdot = vertcat(dtrajectory)

_tv_p = vertcat(tv_param_1)
"""
--------------------------------------------------------------------------
template_model: initial condition and constraints
--------------------------------------------------------------------------
"""
# Initial conditions for Differential States
i_init  = 0.1
theta_init  = pi/10
v_init = 10

x0 = NP.array([i_init, theta_init, v_init])

# Bounds on the states. Use "inf" for unconstrained states
theta_lb = -180;            theta_ub = 180;
i_lb = -10;                 i_ub =  10;

s_lb = NP.array([theta_lb,i_lb])
s_ub = NP.array([theta_ub,i_ub])

# Bounds on the control inputs. Use "inf" for unconstrained inputs
v_lb=0;             v_ub=100;

u_lb=NP.array([v_lb])
u_ub=NP.array([v_ub])
u0 = NP.array([v_init])

# Scaling factors for the states and control inputs. Important if the system is ill-conditioned
x_scaling = NP.array([1.0, 1.0])
z_scaling = NP.array([])
u_scaling = NP.array([1.0])

# Other possibly nonlinear constraints in the form cons(x,u,p) <= cons_ub
# Define the expresion of the constraint (leave it empty if not necessary)
cons = vertcat([])
# Define the lower and upper bounds of the constraint (leave it empty if not necessary)
cons_ub = NP.array([])

# Activate if the nonlinear constraints should be implemented as soft constraints
soft_constraint = 0
# Penalty term to add in the cost function for the constraints (it should be the same size as cons)
penalty_term_cons = NP.array([])
# Maximum violation for the constraints
maximum_violation = NP.array([0])

# Define the terminal constraint (leave it empty if not necessary)
cons_terminal = vertcat([])
# Define the lower and upper bounds of the constraint (leave it empty if not necessary)
cons_terminal_lb = NP.array([])
cons_terminal_ub = NP.array([])

"""
--------------------------------------------------------------------------
template_model: cost function
--------------------------------------------------------------------------
"""
# Define the cost function
# Lagrange term
##-----Distance to vertical position---------------------
lterm =  (v-0.5)**2 +(s-3)**2
# Mayer term
mterm =  0

# Penalty term for the control movements
rterm = 0.0001*NP.array([1.0])

"""
--------------------------------------------------------------------------
template_model: pass information (not necessary to edit)
--------------------------------------------------------------------------
 """
model_dict = {'x':_x,'u': _u, 'rhs':_xdot,'p': _p, 'z':_z,'x0': x0,'x_lb': s_lb,'x_ub': s_ub, 'u0':u0, 'u_lb':u_lb, 'u_ub':u_ub, 'x_scaling':x_scaling, 'u_scaling':u_scaling, 'cons':cons,
"cons_ub": cons_ub, 'cons_terminal':cons_terminal, 'cons_terminal_lb': cons_terminal_lb,'tv_p':_tv_p, 'cons_terminal_ub':cons_terminal_ub, 'soft_constraint': soft_constraint, 'penalty_term_cons': penalty_term_cons, 'maximum_violation': maximum_violation, 'mterm': mterm,'lterm':lterm, 'rterm':rterm}

model = core_do_mpc.model(model_dict)

return model
4flixt commented 4 years ago

do-mpc version 4.0.0 is a complete overhaul of the previous generation. This issue refers to the previous generation and is therefore closed. If the problem persists, please raise a new issue.