How to calculate a multiple correlation with non-negative constraints on the linear model's parameters? From wikipedia:
https://en.wikipedia.org/wiki/Multiple_correlation

In statistics, the coefficient of multiple correlation is a measure of
how well a given variable can be predicted using a linear function of
a set of other variables.

Is there a modification (standard names and methodology, python software) of multiple correlation where linear function is replaced with conic function? I.e., a linear function allows negative coefficients and conic function does not.

One possibility I could think of is to use mathematical programming to resolve this......
 A: If I understand this right, you can estimate a multiple regression model with non-negativity restrictions on the coefficients (in R, this can be done with, for instance,the CRAN package nnls), and then use the R-squared from that fit.  There might well be some similar functions in python.
A: Core Answer
Echoing the answer by Kjetil, you could approach this using non-negative least squares followed by calculating the $R^2$ for the fitted model.
In Python you can use scipy.optimize.nnls.
Example 1
Here is an example usage adapted from the documentation:
import numpy as np
from scipy.optimize import nnls

# Make up some data
m, n = 100, 2
X = np.random.normal(size=m*n).reshape(m, n)
theta = np.arange(2) + 1
y = X @ theta

# Fit model
params, ss_res =  nnls(X, y)

# Compute R^2
r_squared = 1 - ss_res / np.var(y) 

Note that above I did not include an intercept, and as a consequence of that the linear model we trained could have a negative $R^2$ (although not on the idealized data we used in this case).
Example 2
This second example includes an intercept, which is achieved by including a column of ones and an additional parameter.
import numpy as np
from scipy.optimize import nnls

# Make up some data
m, n = 100, 2
X = np.random.normal(size=m*n).reshape(m, n)
X = np.concatenate((X, np.ones(m).reshape(-1, 1)), axis=1)
theta = np.arange(3) + 1
y = X @ theta

# Fit model
params, ss_res =  nnls(X, y)

# Compute R^2
r_squared = 1 - ss_res / np.var(y)

