Skip to content

Soft constraints in AVIs #108

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
darnstrom merged 2 commits into from
Jan 18, 2026
Merged

Soft constraints in AVIs #108

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
e7c6499
Soft constraints for AVIs
darnstrom Jan 17, 2026
93574dd
Bump version
darnstrom Jan 17, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion interfaces/daqp-julia/Project.toml
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
name = "DAQPBase"
uuid = "833476c3-a8c0-4073-9b64-6473509843fe"
authors = ["Daniel Arnström <daniel.arnstrom@gmail.com>"]
version = "0.3.5"
version = "0.3.6"

[deps]
DAQP_jll = "5c51c916-43bf-57fe-9d62-13064ebbf40d"
Expand Down
115 changes: 54 additions & 61 deletions interfaces/daqp-julia/src/avi.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,16 +10,20 @@ end
struct AVIWorkspace
H::Matrix{Float64}
f::Vector{Float64}
A::Matrix{Float64}
At::Matrix{Float64}
bu::Vector{Float64}
bl::Vector{Float64}
sense::Vector{Cint}
norm_factors::Vector{Cdouble}

H_factor::LU{Float64}
rho_soft::Float64

H_factor::LU{Float64, Matrix{Float64}, Vector{Int64}}
x_unc::Vector{Float64}

H1pI::Matrix{Float64}
H2mI::Matrix{Float64}
H2pI::LU{Float64} # LU
H2pI::LU{Float64, Matrix{Float64}, Vector{Int64}}
H2::Matrix{Float64}
H2xpf::Vector{Float64}
x::Vector{Float64}
Expand All @@ -30,13 +34,13 @@ struct AVIWorkspace
settings::AVISettings
end

AVIWorkspace() = AVIWorkspace(zeros(0,0),zeros(0),zeros(0,0),zeros(0),zeros(0),
lu(zeros(0,0)),zeros(0),
AVIWorkspace() = AVIWorkspace(zeros(0,0),zeros(0),zeros(0,0),zeros(0),zeros(0),zeros(Cint,0),zeros(0),
0,lu(zeros(0,0)),zeros(0),
zeros(0,0),zeros(0,0),lu(zeros(0,0)),zeros(0,0), zeros(0),
zeros(0),zeros(0), Model(), zeros(0), AVISettings())


function setup_avi(H,f,A,bu,bl,sense;x0 = zeros(0), daqp_workspace=nothing, settings=AVISettings())
function setup_avi(H,f,A,bu,bl,sense;x0 = zeros(0), rho_soft = 1e6, daqp_workspace=nothing, settings=AVISettings())
n,m = length(f),length(bu)
bl = isempty(bl) ? fill(-1e30,m) : bl
sense = isempty(sense) ? zeros(Cint,m) : sense
Expand All @@ -51,26 +55,52 @@ function setup_avi(H,f,A,bu,bl,sense;x0 = zeros(0), daqp_workspace=nothing, sett
x = isempty(x0) ? fill(NaN,n) : copy(x0)
y,H2xpf = zeros(n),zeros(n)

kkt_buffer = Float64[]
sizehint!(kkt_buffer,4*(n+1)*(n+1))
kkt_buffer = zeros(4*(n+1)*(n+1))
#sizehint!(kkt_buffer,4*(n+1)*(n+1))

daqp_workspace = isnothing(daqp_workspace) ? Model() : daqp_workspace;
DAQPBase.settings(daqp_workspace, Dict(:iter_limit => settings.inner_iter_limit,
:primal_tol => settings.primal_tol,
:dual_tol => settings.dual_tol))
exitflag,tsetup = DAQPBase.setup(daqp_workspace, H1pI, H2xpf, A, bu, bl, sense)

workspace = unsafe_load(daqp_workspace.work);
norm_factors = unsafe_wrap(Vector{Cdouble}, workspace.scaling, m; own = false)


if m > size(A,1)
A = [I(n)[1:m-size(A,1),:];A] # To handle simple bounds
#A = [I(n)[1:m-size(A,1),:];A] # To handle simple bounds
At = [I(n)[:,1:m-size(A,1)] A'] # To handle simple bounds
else
At = A'[:,:]
end

return exitflag, AVIWorkspace(H,f,A,bu,bl,
H_factor,zeros(n),
return exitflag, AVIWorkspace(H,f,At,bu,bl,sense,norm_factors,
rho_soft,H_factor,zeros(n),
H1pI,H2mI,H2pI,H2,H2xpf,
x,y,daqp_workspace,kkt_buffer,
settings)
end

function _is_optimal(λ,xt,λfull,ws,ASu,ASl)
nu = length(ASu)
ϵp, ϵd = ws.settings.primal_tol, ws.settings.dual_tol
@inbounds for i in 1:nu
λ[i] < -ϵd && return false
end
@inbounds for i in nu+1:length(λ)
λ[i] > ϵd && return false
end
@inbounds for i in 1:length(ws.bu)
if λfull[i] == 0.0 # Only test inactive constraints
Axi = dot(view(ws.At,:,i),xt)
Axi-ws.bu[i] > ϵp && return false
Axi-ws.bl[i] < -ϵp && return false
end
end
return true
end

function solve(ws::AVIWorkspace)
exitflag = 0
λstar, ASstar = zeros(0), Int[]
Expand All @@ -91,17 +121,12 @@ function solve(ws::AVIWorkspace)
if info.iterations == 1
nkkt += 1
ASu,ASl = _get_AS(info.λ,ϵd)
#xt,λ,AS = _solve_kkt_schur(ws,ASu,ASl)
xt,λ,AS = _solve_kkt(ws,ASu,ASl)
nu = length(ASu)
if all(λ[i] > -ϵd for i in 1:nu) && all(λ[i] < ϵd for i in nu+1:length(λ))
Ax = ws.A*xt
if all(Ax-ws.bu .<= ϵp) && all(Ax-ws.bl .>=-ϵp)
ws.x .= xt
λstar,ASstar = copy(λ),AS
exitflag, outer_iter = 1,k
break
end
if _is_optimal(λ,xt,info.λ,ws,ASu,ASl)
ws.x .= xt
λstar,ASstar = copy(λ),AS
exitflag, outer_iter = 1,k
break
end
axpby!(β,xt,1.0-β,ws.x)
else
Expand All @@ -127,53 +152,21 @@ function _get_AS(λ,dual_tol)
return ASu,ASl
end

# TODO kkt_schur is faster in isolation, but leads to slower total execution
#function _solve_kkt_schur(ws::AVIWorkspace,ASu,ASl)
# # Prepare some helpers
# AS = [ASu;ASl]
# n,nAS = length(ws.f),length(AS)
#
# # Juggle some buffers
# resize!(ws.kkt_buffer,nAS*(nAS+n+1)+n) # Is this necessary?
# invH_At = reshape(view(ws.kkt_buffer,1:nAS*n),(n,nAS))
# offset = nAS*n
# S = reshape(view(ws.kkt_buffer,offset+1:offset+nAS^2),(nAS,nAS))
# offset +=nAS^2
# λ = view(ws.kkt_buffer,offset+1:offset+nAS)
# offset +=nAS
# x = view(ws.kkt_buffer,offset+1:offset+n)
#
# AAS = ws.A[AS,:]
#
# # Schur complement
# @inbounds for (i,id) in enumerate(AS)
# @views invH_At[:,i] .= ws.A[id,:]
# end
# ldiv!(ws.H_factor,invH_At)
# mul!(S, AAS, invH_At)
# S_fact = lu!(S)
#
# # Solve for λ
# @views λ .= AAS * ws.x_unc
# @views λ[1:length(ASu)] .-= ws.bu[ASu]
# @views λ[length(ASu)+1:end] .-= ws.bl[ASl]
# ldiv!(S_fact,λ)
#
# # Solve for x
# x .= ws.x_unc .- invH_At * λ
#
# return x,λ,AS
#end
function _solve_kkt(ws,ASu,ASl)
AS = [ASu;ASl]
AS = Int[ASu;ASl]
n = length(ws.f)
nkkt = n+length(AS)
resize!(ws.kkt_buffer,nkkt^2)
K = reshape(view(ws.kkt_buffer,1:nkkt^2),(nkkt,nkkt))
@views K[1:n,1:n] = ws.H
@views K[n+1:end,1:n] = ws.A[AS,:]
@views K[1:n,n+1:end] = ws.A[AS,:]'
@views K[n+1:end,1:n] = ws.At[:,AS]'
@views K[1:n,n+1:end] = ws.At[:,AS]
K[n+1:end,n+1:end] .= 0
for (i,id) in enumerate(AS)
if ws.sense[id] & DAQPBase.SOFT != 0
K[n+i,n+i] = -1/((ws.norm_factors[id]^2)*ws.rho_soft)
end
end
z = K\[-ws.f;ws.bu[ASu];ws.bl[ASl]]
x = @view z[1:n]
λ = @view z[n+1:end]
Expand Down
Loading