search

JuliaMPC / NLOptControl.jl

nonlinear control optimization tool
Other
109 stars 26 forks source link

Simple Optimal control problem #32

Closed azev77 closed 3 years ago

azev77 commented 3 years ago

@huckl3b3rry87 I'm able to get NLOptControl.jl to work now. Thank you!

I'm having trouble solving the following seemingly simple problem: image

using NLOptControl
n=define(
    numStates=1,    # b(t) wealth
    numControls=1,  # c(t) consumption
    X0=[100.0],     # b(0)= 100.0
    XF=[0.0]        # b(10)= 0.0 die w/o assets
    )
states!(n,[:b];descriptions=["b(t)"])
controls!(n,[:c];descriptions=["c(t)"])
dx=[:( (0.05) * b[j] - c[j] )]
dynamics!(n,dx)
configure!(n;(:Nck=>[50]),(:finalTimeDV=>false),(:tf=>10.0));
obj=integrate!(n, :( -1.0 * exp(-.01*j) * log(c[j]) ) )
@NLobjective(n.ocp.mdl,Min,obj)
optimize!(n);
allPlots(n)

#Error message

┌ Warning: Ipopt finished with status Invalid_Number_Detected
└ @ Ipopt ~/.julia/packages/Ipopt/bYzBL/src/MPB_wrapper.jl:178
┌ Warning: Not solved to optimality, status: Error
└ @ JuMP ~/.julia/packages/JuMP/I7whV/src/nlp.jl:1283
┌ Warning: The solution is not Optimal 
│ 
│                 Setting: n.f.mpc.simFailed = [true, n.r.ocp.status] 
└ @ NLOptControl ~/.julia/packages/NLOptControl/i4KEb/src/utils.jl:406

julia> 

# I tried to change the configure options w/o luck
#configure!(n;(:integrationScheme=>:trapezoidal),(:finalTimeDV=>false),(:tf=>10.0));
#configure!(n;(:finalTimeDV=>false),(:tf=>1.0));
#configure!(n;(:tf=>1.0));

Do you know why it's not working?

huckl3b3rry87 commented 3 years ago

What is the error

On Tue, Oct 13, 2020, 3:52 PM azev77 notifications@github.com wrote:

@huckl3b3rry87 https://github.com/huckl3b3rry87 I'm able to get NLOptControl.jl to work now. Thank you!

I'm having trouble solving the following seemingly simple problem: [image: image] https://user-images.githubusercontent.com/7883904/95923957-9dd8b000-0d84-11eb-98ba-f75dbddde9c0.png

using NLOptControl n=define( numStates=1, # b(t) wealth numControls=1, # c(t) consumption X0=[100.0], # b(0)= 100.0 XF=[0.0] # b(10)= 0.0 die w/o assets )states!(n,[:b];descriptions=["b(t)"])controls!(n,[:c];descriptions=["c(t)"]) dx=[:( (0.05) b[j] - c[j] )]dynamics!(n,dx)configure!(n;(:Nck=>[100]),(:finalTimeDV=>false),(:tf=>10.0));#configure!(n;(:integrationScheme=>:trapezoidal),(:finalTimeDV=>false),(:tf=>10.0));#configure!(n;(:finalTimeDV=>false),(:tf=>10.0));#configure!(n;(:tf=>10.0)); obj=integrate!(n, :( -1.0 exp(-.01j) log(c[j]) ) )@NLobjective(n.ocp.mdl,Min,obj)optimize!(n);allPlots(n)

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/JuliaMPC/NLOptControl.jl/issues/32, or unsubscribe https://github.com/notifications/unsubscribe-auth/AB6TAWWP2AI65RZUUALB3ADSKTK4PANCNFSM4SPZ53JQ .

azev77 commented 3 years ago

@huckl3b3rry87 the error message is

┌ Warning: Ipopt finished with status Invalid_Number_Detected
└ @ Ipopt ~/.julia/packages/Ipopt/bYzBL/src/MPB_wrapper.jl:178
┌ Warning: Not solved to optimality, status: Error
└ @ JuMP ~/.julia/packages/JuMP/I7whV/src/nlp.jl:1283
┌ Warning: The solution is not Optimal 
│ 
│                 Setting: n.f.mpc.simFailed = [true, n.r.ocp.status] 
└ @ NLOptControl ~/.julia/packages/NLOptControl/i4KEb/src/utils.jl:406
huckl3b3rry87 commented 3 years ago

This is an issue with your formulation. It is likely due to a divide by zero error or something . I would suggest to think about when the formulation is undefined, for instance, try constraining c to be greater than 0.1 or something.

On Tue, Oct 13, 2020, 6:47 PM azev77 notifications@github.com wrote:

@huckl3b3rry87 https://github.com/huckl3b3rry87 the error message is

┌ Warning: Ipopt finished with status Invalid_Number_Detected

└ @ Ipopt ~/.julia/packages/Ipopt/bYzBL/src/MPB_wrapper.jl:178

┌ Warning: Not solved to optimality, status: Error

└ @ JuMP ~/.julia/packages/JuMP/I7whV/src/nlp.jl:1283

┌ Warning: The solution is not Optimal

│

│ Setting: n.f.mpc.simFailed = [true, n.r.ocp.status]

└ @ NLOptControl ~/.julia/packages/NLOptControl/i4KEb/src/utils.jl:406

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/JuliaMPC/NLOptControl.jl/issues/32#issuecomment-708103616, or unsubscribe https://github.com/notifications/unsubscribe-auth/AB6TAWWRSEA4ONRY7AQHRZLSKT7KJANCNFSM4SPZ53JQ .

azev77 commented 3 years ago

@huckl3b3rry87 I tried constraining c > 0.1. No luck unfortunately...

huckl3b3rry87 commented 3 years ago

Which version of NLOptControl are you using? For the newest version, I am not sure about pseudospectral methods, try trapezoidal.

And please provide full input/output

On Tue, Oct 13, 2020, 9:04 PM azev77 notifications@github.com wrote:

@huckl3b3rry87 https://github.com/huckl3b3rry87 I tried constraining c > 0.1. No luck unfortunately...

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/JuliaMPC/NLOptControl.jl/issues/32#issuecomment-708142631, or unsubscribe https://github.com/notifications/unsubscribe-auth/AB6TAWV6KE4MOCFBCXVXLKDSKUPODANCNFSM4SPZ53JQ .

huckl3b3rry87 commented 3 years ago

Try constraining both the state and the control with upper and lower bounds.

On Tue, Oct 13, 2020, 9:12 PM Huckleberry Febbo febbo@umich.edu wrote:

Which version of NLOptControl are you using? For the newest version, I am not sure about pseudospectral methods, try trapezoidal.

And please provide full input/output

On Tue, Oct 13, 2020, 9:04 PM azev77 notifications@github.com wrote:

@huckl3b3rry87 https://github.com/huckl3b3rry87 I tried constraining c > 0.1. No luck unfortunately...

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/JuliaMPC/NLOptControl.jl/issues/32#issuecomment-708142631, or unsubscribe https://github.com/notifications/unsubscribe-auth/AB6TAWV6KE4MOCFBCXVXLKDSKUPODANCNFSM4SPZ53JQ .

azev77 commented 3 years ago

The code I tried:

using NLOptControl
n=define(
    numStates=1,    # b(t) wealth
    numControls=1,  # c(t) cons
    X0=[10.0],     # b(0)= 100.0
    XF=[0.0],      # b(tf)= 0.0 die w/o assets
    CL=[1.1],CU=[100.],      # make sure cons is >0
    XL=[-50.0],XU=[50.0]
    )
states!(n,[:b];descriptions=["b(t)"])
controls!(n,[:c];descriptions=["c(t)"])
dx=[:( (0.05) * b[j] - c[j] )]
dynamics!(n,dx)
# configure!(n;(:Nck=>[1000]),(:finalTimeDV=>false),(:tf=>10.0));
configure!(n;(:integrationScheme=>:trapezoidal),(:finalTimeDV=>false),(:tf=>10.0));
#configure!(n;(:finalTimeDV=>false),(:tf=>1.0));
#configure!(n;(:tf=>1.0));
obj=integrate!(n, :( -1.0 * exp(-.01*j) * log(c[j]) ) )
@NLobjective(n.ocp.mdl, Min, obj)
optimize!(n)

The Error Message:

 Warning: Ipopt finished with status Invalid_Number_Detected
└ @ Ipopt ~/.julia/packages/Ipopt/bYzBL/src/MPB_wrapper.jl:178
┌ Warning: Not solved to optimality, status: Error
└ @ JuMP ~/.julia/packages/JuMP/I7whV/src/nlp.jl:1283
┌ Warning: The solution is not Optimal 
│ 
│                 Setting: n.f.mpc.simFailed = [true, n.r.ocp.status] 
└ @ NLOptControl ~/.julia/packages/NLOptControl/i4KEb/src/utils.jl:406
huckl3b3rry87 commented 3 years ago

Try playing with the bounds on state and control, reduce them a bunch

On Tue, Oct 13, 2020, 9:36 PM azev77 notifications@github.com wrote:

The code I tried:

using NLOptControl

n=define(

numStates=1,    # b(t) wealth

numControls=1,  # c(t) cons

X0=[10.0],     # b(0)= 100.0

XF=[0.0],      # b(tf)= 0.0 die w/o assets

CL=[1.1],CU=[100.],      # make sure cons is >0

XL=[-50.0],XU=[50.0]

)

states!(n,[:b];descriptions=["b(t)"]) controls!(n,[:c];descriptions=["c(t)"])

dx=[:( (0.05) * b[j] - c[j] )] dynamics!(n,dx)

configure!(n;(:Nck=>[1000]),(:finalTimeDV=>false),(:tf=>10.0));

configure!(n;(:integrationScheme=>:trapezoidal),(:finalTimeDV=>false),(:tf=>10.0));

configure!(n;(:finalTimeDV=>false),(:tf=>1.0));

configure!(n;(:tf=>1.0));

obj=integrate!(n, :( -1.0 exp(-.01j) * log(c[j]) ) ) @NLobjective(n.ocp.mdl, Min, obj) optimize!(n)

The Error Message:

Warning: Ipopt finished with status Invalid_Number_Detected

└ @ Ipopt ~/.julia/packages/Ipopt/bYzBL/src/MPB_wrapper.jl:178

┌ Warning: Not solved to optimality, status: Error

└ @ JuMP ~/.julia/packages/JuMP/I7whV/src/nlp.jl:1283

┌ Warning: The solution is not Optimal

│

│ Setting: n.f.mpc.simFailed = [true, n.r.ocp.status]

└ @ NLOptControl ~/.julia/packages/NLOptControl/i4KEb/src/utils.jl:406

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/JuliaMPC/NLOptControl.jl/issues/32#issuecomment-708151009, or unsubscribe https://github.com/notifications/unsubscribe-auth/AB6TAWSGC6CJEYWS7PCYKDLSKUTGHANCNFSM4SPZ53JQ .

azev77 commented 3 years ago

Thanks! I was ready to give up, but now it works!

using NLOptControl
n=define(
    numStates=1, numControls=1,    
    X0=[5.0], XF=[0.0],     # b(0)= 100.0
    CL=[0.0], CU=[50.0],      # make sure cons is >0
    XL=[-50.0],XU=[50.0]
    )
states!(n,[:b];descriptions=["b(t)"])
controls!(n,[:c];descriptions=["c(t)"])
dx=[:( (0.01) * b[j] + 10.0 - c[j] )]
dynamics!(n,dx)
configure!(n;(:integrationScheme=>:trapezoidal),(:finalTimeDV=>false),(:tf=>10.0))
obj=integrate!(n, :( exp(-.04 * j) * log(c[j]) ) );
@NLobjective(n.ocp.mdl, Max, obj)
optimize!(n)
allPlots(n)
#
##########
huckl3b3rry87 commented 3 years ago

Awesome 👍

On Tue, Oct 13, 2020, 10:14 PM azev77 notifications@github.com wrote:

Thanks! I was ready to give up, but now it works!

using NLOptControl n=define( numStates=1, numControls=1, X0=[5.0], XF=[0.0], # b(0)= 100.0 CL=[0.0], CU=[50.0], # make sure cons is >0 XL=[-50.0],XU=[50.0] )states!(n,[:b];descriptions=["b(t)"])controls!(n,[:c];descriptions=["c(t)"]) dx=[:( (0.01) b[j] + 10.0 - c[j] )]dynamics!(n,dx)configure!(n;(:integrationScheme=>:trapezoidal),(:finalTimeDV=>false),(:tf=>10.0)) obj=integrate!(n, :( exp(-.04 j) * log(c[j]) ) );@NLobjective(n.ocp.mdl, Max, obj)optimize!(n)allPlots(n)###########

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/JuliaMPC/NLOptControl.jl/issues/32#issuecomment-708161217, or unsubscribe https://github.com/notifications/unsubscribe-auth/AB6TAWQBXXI7WWPFKTE3NXTSKUXTBANCNFSM4SPZ53JQ .

azev77 commented 3 years ago

@huckl3b3rry87

  1. is it possible to set an inequality constraint so that a control is less than or equal to a state var? u1[j] <= 0.25 * x2[j]
  2. is it possible to solve an infinite horizon model tf=infinity?

I can't find how in the docs?

huckl3b3rry87 commented 3 years ago

  1. Yes, use the @NLConstraint JuMP macro

  2. No. Finite time only

On Fri, Oct 16, 2020, 10:07 PM azev77 notifications@github.com wrote:

@huckl3b3rry87 https://github.com/huckl3b3rry87

  1. is it possible to set an inequality constraint so that a control is less than or equal to a state var? u1[j] <= 0.25 * x2[j]
  2. is it possible to solve an infinite horizon model tf=infinity?

I can't find how in the docs?

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/JuliaMPC/NLOptControl.jl/issues/32#issuecomment-710726795, or unsubscribe https://github.com/notifications/unsubscribe-auth/AB6TAWWFPCFLUA5ZJKO6USDSLD34HANCNFSM4SPZ53JQ .

azev77 commented 3 years ago

I can't find examples using @NLConstraint in NLOptControl:

using NLOptControl;
n=define(
    numStates=2, numControls=2,
    X0=[5.0,0.0], XF=[0.01,0.0],
    XL=[0.0,NaN], XU=[100.0,NaN],
    CL=[0.01,-7.0], CU=[50.0,NaN]
    )
states!(n,[:k,:b];descriptions=["k(t)","b(t)"])
controls!(n,[:c,:pmt];descriptions=["c(t)","pmt(t)"])
dx=[
    :( (k[j]^(.75)) - 0.3 * k[j] - c[j] - pmt[j] ),
    :( 0.01 * b[j] - pmt[j] )
    ]
dynamics!(n,dx)
configure!(n;
    (:integrationScheme=>:trapezoidal),
    (:finalTimeDV=>false),(:tf=>100.0))
#
@NLconstraint(n.ocp.NLcon, :( b[j]  <= 0.5 * k[j]) )    # ERROR 
#
obj=integrate!(n, :( exp(-.03 * j) * log(c[j]) ) );
@NLobjective(n.ocp.mdl, Max, obj)
optimize!(n)
allPlots(n)

What is the correct syntax for: @NLconstraint(n.ocp.NLcon, :( b[j] <= 0.5 * k[j]) ) ?

huckl3b3rry87 commented 3 years ago

Look through this package https://github.com/JuliaMPC/MichiganAutonomousVehicles.jl

for some more complicated examples. For instance https://github.com/JuliaMPC/MichiganAutonomousVehicles.jl/blob/d4764c100e5e9a5d309e3e236b2fc9fadd0f24cb/src/AutonomousControl.jl#L396

Let me know if you figure it out

On Sat, Oct 17, 2020, 1:06 PM azev77 notifications@github.com wrote:

I can't find examples using @NLConstraint in NLOptControl:

using NLOptControl; n=define( numStates=2, numControls=2, X0=[5.0,0.0], XF=[0.01,0.0], XL=[0.0,NaN], XU=[100.0,NaN], CL=[0.01,-7.0], CU=[50.0,NaN] )states!(n,[:k,:b];descriptions=["k(t)","b(t)"])controls!(n,[:c,:pmt];descriptions=["c(t)","pmt(t)"]) dx=[ :( (k[j]^(.75)) - 0.3 k[j] - c[j] - pmt[j] ), :( 0.01 b[j] - pmt[j] ) ]dynamics!(n,dx)configure!(n; (:integrationScheme=>:trapezoidal), (:finalTimeDV=>false),(:tf=>100.0))#@NLconstraint(n.ocp.NLcon, :( b[j] <= 0.5 k[j]) ) # ERROR # obj=integrate!(n, :( exp(-.03 j) * log(c[j]) ) );@NLobjective(n.ocp.mdl, Max, obj)optimize!(n)allPlots(n)

What is the correct syntax for: @NLconstraint(n.ocp.NLcon, :( b[j] <= 0.5

  • k[j]) )

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/JuliaMPC/NLOptControl.jl/issues/32#issuecomment-711048338, or unsubscribe https://github.com/notifications/unsubscribe-auth/AB6TAWSU2DOQC4JYPHGB6ALSLHFH7ANCNFSM4SPZ53JQ .